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Abstract 

We revisit Lyman-a bounds on the dark matter mass in A Warm Dark 
Matter (AWDM) models, and derive new bounds in the case of mixed Cold plus 
Warm models (AC WDM), using a set up which is a good approximation for 
several theoretically well-motivated dark matter models. We combine WMAP5 
results with two different Lyman-a data sets, including observations from the 
Sloan Digital Sky Survey. We pay a special attention to systematics, test 
various possible sources of error, and compare the results of different statistical 
approaches. Expressed in terms of the mass of a non-resonantly produced 
sterile neutrino, our bounds read ?nNRp > 8 keV (frequentist 99.7% confidence 
limit) or ttinrp > 12.1 keV (Bayesian 95% credible interval) in the pure AWDM 
limit. For the mixed model, we obtain limits on the mass as a function of 
the warm dark matter fraction .Fwdm- Within the mass range studied here 
(5 keV < ttinrp < oo), we find that any mass value is allowed when Fwdm < 0.6 
(frequentist 99.7% confidence limit); similarly, the Bayesian joint probability 
on (Fwdm, 1/'71nrp) allows any value of the mass at the 95% confidence level, 
provided that Fwdm < 0.35. 
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1 Introduction 



The nature of Dark Matter (DM) is one of the most intriguing questions of particle 
astrophysics. Its resolution would have a profound impact on the development of 
particle physics beyond its Standard Model (SM). Although the possibility of having 
massive compact halo objects (MACHOs) as a dominant form of DM is still under 
debates (see recent discussion in [1] and references therein), it is widely believed that 
DM is made of non-baryonic particles. Yet the SM of elementary particles does not 
contain a viable DM particle candidate - a massive, neutral and stable (or at least 
cosmologically long-lived) particle. Active neutrinos, which are both neutral and 
stable, form structures in a top-down fashion [2, 3, 4, 5, 6], and thus cannot produce 
the observed large scale structure. Therefore, the DM particle hypothesis implies an 
extension of the SM. Thus, constraining the properties of DM helps to distinguish 
between various DM candidates and may help to differentiate among different models 
Beyond the Standard Model (BSM). 

What is known about the properties of DM particles? 

DM particle candidates may have very different masses (for reviews of DM can- 
didates see e.g. [7, 8, 9, 10]): massive gravitons with mass ~ 10""^^ eV [11], ax- 
ions with mass ~ 10~^ eV [12, 13, 14], sterile neutrinos with mass in the keV 
range [15], supersymmetric (SUSY) particles (gravitinos [16], neutralinos [17, 18], 
axinos [19, 20] with mass ranging from a few eV to hundreds GeV), supersymmet- 
ric Q-balls [21], WIMPZILLAs with mass ~ 10^^ GeV [22, 23], and many others. 
Thus, the mass of DM particles becomes an important characteristics which may 
help to distinguish between various DM candidates and, more importantly, may help 
to differentiate among different BSMs. A quite robust and model-independent lower 
bound on the mass of spin-one-half DM particles was suggested in [2 1] . The idea 
was based on the fact that for any fermionic DM the average phase-space density 
(in a given DM-dominated, gravitationally bound object) cannot exceed the phase- 
space density of the degenerate Fermi gas. This argument, applied to the most 
DM-dominated dwarf Spheroidal (dSph) satellites of the Milky Way leads to the 
bound > 0.41 keV [25]. For particular DM models (with known primordial ve- 
locity dispersion) and under certain assumptions about the evolution of the system 
during structure formation, this limit can be strengthened. This idea was developed 
in a number of works [26, 27, 28, 29, 30, 31, 32, 33, 25, 34] (for recent results and a 
critical discussion see [25]). 

For any DM candidate there should exist a production mechanism in the early 
Universe. Although it is possible that the DM is produced entirely through interac- 
tions with non-SM particles (e.g. from the infiaton decay) and is inert with respect to 
all SM interactions, many popular DM candidates are produced via interactions with 
the SM sector. One possible class of DM candidates are weakly interacting massive 
particles - WIMPs [17]. They have interactions with SM particles of roughly elec- 
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troweak strength. The WIMPs are completely stable particles, as their interaction 
strength would otherwise imply a lifetime incompatible with a sizable relic density 
today. Therefore, they are produced in the early Universe via annihilation into SM 
particles with roughly weak interaction cross-section (see e.g. [35]). To have a correct 
abundance, these particles should have a mass in the GeV to TeV range. Such a mass 
and interaction strength means that they decouple from the primeval plasma deeply 
in the radiation dominated epoch, while being non-relativistic. Such particles have a 
Boltzmann velocity distribution function at the moment of freeze-out and are called 
"cold DM" (CDM). In the CDM scenario, structure formation goes in the bottom-up 
fashion: the smallest structures (with supposed masses as small as that of the Earth) 
form first, and then combined into larger structures. 

Another wide class of DM candidates can be generically called superweakly inter- 
acting (see e.g. [3G]), meaning that their cross-section of interaction with SM particles 
is much smaller than the weak one. There are many examples of super- WIMP DM 
models: sterile neutrinos [15], gravitinos in theories with broken R-parity [37, 38], 
light volume moduli [39] or Majorons [10]. Many of these candidates (e.g. sterile 
neutrino, gravitino, Majoron) possess a two-body decay channel: DM — > 7 -|- z/, 7 + 7. 
These scenarios are very interesting, since apart from laboratory detection such parti- 
cles can be searched in space. Looking for a monochromatic decay line in the spectra 
of DM-dominated objects provides a way of indirect detection of DM, and could help 
to constrain its interaction strength with SM particles. The astrophysical search for 
decaying DM is in fact more promising than for stable (annihilating) DM. Indeed, 
the decay signal is proportional to the column density f poM{r)dr, i.e. the density 
integrated along the line of sight, and not the squared density / PdmI^)*^^ (^s it is 
the case for annihilating DM). As a consequence, (i) a. vast variety of astrophysical 
objects of different nature would produce roughly the same decay signal [41, 42]; 
(a) this gives more freedom for choosing the observational targets, and avoid com- 
plicated astrophysical backgrounds (e.g. one does not need to look at the Galactic 
center, expecting a comparable signal from dark outskirts of galaxies, clusters and 
dark dSph's); (Hi) if a candidate line is identified, its surface brightness profile may 
be measured (since it does not decay quickly away from the centers of objects), in 
order to distinguish it from astrophysical lines (which usually decay in outskirts) 
and to compare with the signal from several other objects. For all these reasons, 
astrophysical search for decaying DM can almost be viewed as another type of direct 
detection. A search for DM decay signals was conducted both in the keV - MeV 
range [43, 44, 41, 45, 46, 47] and GeV range [42]. 

The small strength of interaction of light super- WIMP particles usually implies 
that: (i) they were never in thermal equilibrium in the early Universe; and (ii) they 
were produced deep in the radiation dominated epoch, while still being relativistic 
(unlike CDM). The latter property means that the density perturbation of these 
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particles are suppressed for scales below their free-streaming length, which affects 
structure formation on those scales. If the maximum free-streaming length roughly 
coincides with galaxy scales, the particles are called Warm Dark Matter (WDM), see 
e.g. [48]. 

Currently, the best way to distinguish between WDM and CDM models is the 
analysis of Lyman-a (Ly-a) forest data (for an introduction see e.g. [49, 50, 51]). 
This method requires particular care. Part of the physics entering into the Lj-a 
analysis is not fully understood (for a recent review see [52]), and can be significantly 
influenced by DM particles (see e.g. [53, 54, 55, 56, 57, 58, 59]). The results can 
also be affected by approximations related to computational difficulties. Indeed, 
at the redshifts probed by Ly-a data, structures are already undergoing non-linear 
evolution. In order to relate the measured flux spectrum with the parameters of 
each cosmological model, one would have to perform a prohibitively large number 
of numerical simulations. Therefore, various simplifying approximations have to be 
realized. Different approaches are presented in [60, 61, 62]. 

Previous works analyzing Lj-a constrains on WDM [63, 65, 64, 66, 67, 68] as- 
sumed the simplest AWDM model with a cut-off in the linear power spectrum (PS) 
of matter density fluctuations. Such a PS arises when WDM particles are ther- 
mal relics. However, in many super-weakly interacting DM models, the PS has a 
complicated non-universal form due to the non-thermal primordial velocity distri- 
bution (see e.g. [69, 70, 71, 72]). For example, in gravitino models with broken 
R-parity, gravitino production occurs in two stages: thermally at high temperatures 
(see e.g. [73, 71, 75]), and then non-thermally through the late decay of the next-to- 
lightest supersymmetric particle (see e.g. [76, 69]). Therefore the primordial velocity 
distribution results from the superimposition of a colder (resonant) and a warmer 
(thermal) component, and the PS is characterized by a step and a plateau on small 
scales (rather than a sharp cut-off). 

Another interesting super- WIMP WDM model with a non-trivial PS is the sterile 
neutrino DM of the i^MSM- an extension of the SM by three sterile (right-handed, 
gauge singlet) neutrinos [77, 78]. Having masses of all new particles below the electro- 
weak scale, this model is able to explain simultaneously three BSM phenomena: 
(i) the transition between the neutrino of different flavours {neutrino oscillations); 
(a) the production of a non-zero baryon number in the early Universe; and (Hi) the 
existence of DM candidate with correct relic density.^ Sterile neutrino DM in the 
z/MSM is produced via its mixing with active neutrinos in the early Universe [15, 

^The free-streaming length is the equivalent of the Jeans length for a coUisionless fluid. It is 
given by the ratio Xps ^ where (v) is the velocity dispersion of the particles. 

^This model also eases the hierarchy problem, as the z/MSM can be thought of as a consistent 
quantum field theory, valid all the way to the Planck scale. The hierarchy problem is then shifted 
into the realm of quantum gravity, broadly understood as a fundamental theory, superseeding at 
the Planck scale [79]. 
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80, 81, 82, 83]. In the absence of lepton asymmetry, the production rate is strongly 
suppressed at temperatures above few hundreds of MeV and peaks roughly at Tpeak ~ 

130 (y^y) MeV [15]. The resulting momentum distribution of sterile neutrino DM 
is approximately that of a rescaled Fermi-Dirac [15, 84], with its temperature equal 
to that of active neutrinos (see (1) below). We will call this production mechanism 
NRP - non-resonant production. Notice that this production channel is always 
present (although it may be subdominant, see below). 

The sterile neutrino DM production changes in the presence of a lepton asym- 
metry [80]. In this case the dispersion relations for active and sterile neutrinos may 
cross each other. This results in an effective transfer of an excess of active neutrinos 
(or antineutrinos) to the population of DM sterile neutrinos. The lepton asymmetry 
necessary for this resonant production (RP) is generated via the decay of heavier 
sterile neutrinos [85]. In the RP scenario, for each mass and lepton asymmetry, the 
velocity dispersion has a colder (resonant) component and a warmer (non-resonantly 
produced) tail. For a recent review see e.g. [86]. 

Sterile neutrinos having colder and warmer velocity components can also be pro- 
duced through the decay of a gauge-singlet scalar, followed by a subsequent non- 
resonant production [87, 88, 89, 90, 91] (see also [92, 93] for other mechanisms of 
production of sterile neutrino). The decay of this scalar field should occur at tem- 
peratures > 100 GeV and the produced sterile neutrinos do not share subsequent 
entropy releases [87, 89]. As a result the characteristic momentum of sterile neutri- 
nos, produced in that way, is several times colder than the average momentum of the 
NRP component (see e.g. [87]). Notice that in this case the production is determined 
not only by the mixing between active and sterile neutrinos, but also by the coupling 
of the sterile neutrino with the scalar. 

We see that in several well-motivated models, the primordial momentum distribu- 
tion of DM particles contains a warm (Fermi-Dirac-like component) and a colder one. 
Therefore, in this work, apart from updating existing bounds on AWDM models, we 
consider Ly-a bounds in the ACWDM scenario, including a mixture of a CDM and 
WDM. This model is characterized by two independent parameters beyond the usual 
dark matter fraction ^dm^ the velocity dispersion (or free-streaming length) of the 
WDM component, and the fraction of WDM density. We will overview and provide 
a critical analysis of uncertainties in the Ly-a method for this model, and derive 
constraints on the above parameters. Finally, we will discuss the consequences of 
these results for some physically relevant DM models of this kind. 

The first qualitative analysis of Ly-a bounds in the ACWDM scenario was per- 
formed in [91]. In this work, the authors computed the linear ACWDM power 
spectra projected along the line of sight. The PS suppression (as compared to the 
pure ACDM case) at some fiducial scale was compared with that quoted in [66] for 
several WDM masses. A ACWDM model was considered as ruled out at some con- 
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fidence level if it produced the same suppression in the ID projected linear PS as a 
pure AWDM model ruled out at the same level in [66]. 

This qualitative method should be improved in two directions. First of all, the 
Ly-a method provides a measurement of the ID non-linear flux spectrum. The non- 
linear evolution "mixes" different scales. Therefore, two linear PS passing through 
the same point at a given wave number but having a different behavior otherwise 
would fit Lj-a data differently. Secondly, the effects of CWDM suppression on the 
matter PS can be somewhat compensated by the change of cosmological parameters 
(DM abundances, erg, etc.). Therefore, a proper Ly-a analysis should include a 
simultaneous fitting of all cosmological plus CWDM parameters to the Ly-a (and 
possibly other cosmological) data. Such an analysis is performed in this work. 

The paper is organized as follows. In section 2, we describe the effect of a mixture 
of cold and warm dark matter on the linear matter power spectrum, and present some 
analytic approximations. In section 3, we introduce the two Ly-a data sets used in 
this analysis: the Viel, Haenelt & Springel (VHS) data compiled in Ref. [95], and 
the Sloan Digital Sky Survey (SDSS) data complied in Ref. [96, 61]. In section 4, 
we discuss the issue of including thermal velocities in hydrodynamical simulations 
in presence of warm dark matter. The reader only interested in our results could go 
directly to section 5, where we present the bounds based on the conservative VHS 
data set, and to section 6, which contains our analysis for the more constraining 
SDSS data. In section 7, we describe the systematic errors included in our analysis. 
Finally, in section 8, we discuss some implications of our results. 



2 Linear power spectrum in CWDM models 

The growth of matter density perturbations in the Universe is conveniently described 
in terms of a power spectrum P{k) = (1^1^) (where pu is the matter density, and 
< ■ > denotes the ensemble average). The time evolution of this power spectrum 
from its primordial form (taken in this paper to be a power law with index rig) to 
a redshift of interest depends on cosmological parameters (fraction of dark fioM and 
baryonic fig matter, Q\, Hubble constant, etc.) 

Pure warm DM models contain another parameter which affects the evolution 
of density perturbations: the initial DM velocity dispersion. In this paper, for def- 
initeness, we assume that the WDM component has a mass mNRp and a primordial 
phase-space distribution f{p) of the form 

= r /'^^ I 1 • 

exp(p/T^) + 1 

This form is motivated by the scenario of non-resonantly produced (NRP) sterile 
neutrinos [15, 84]. The temperature Tj, = (4/ll)^/^T^ is the temperature of active 
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neutrinos in the instantaneous decoupling approximation, while the normalization 
parameter x -C 1 is determined by the mass mNRp and the requirement of a correct 
DM abundance'^. Another well- motivated WDM model is that of thermal relics 
(TR). This is a generic name for particles that were in thermal equilibrium in the 
early Universe and decoupled while being relativistic (unlike sterile neutrinos which 
never equilibrated). In the TR case, the phase-space distribution is of the Fermi- 
Dirac type: namely, similar to equation (1) although x is equal to one, and the 
temperature Ttr is not equal to but can be deduced from the mass of thermal 
relics ttItr and the DM abundance. At the level of linear perturbations, the NRP 
and TR models are formally equivalent under the identification = x^^^T^ ^iid 
mass = X^^^"^nrp [! ' '"]. 

To separate the influence of the velocity dispersion on the evolution of density 
perturbations, it is convenient to construct a transfer function - the square root of 
the ratio of the matter power spectrum of a given AWDM model to that of a pure 
ACDM with the same choice of cosmological parameters: 

1/2 

T{k) 



f -Pawdm a 



This transfer function turns out to be rather insensitive to the values of cosmological 
parameters, and is well fitted by 

T{k) = (1 + {k/kof , (3) 

where v = 1.12, and the scale k^ (related to the free-steaming scale) is given as 
a function of the mass and cosmological parameters (see e.g. [G^)]). This function 
describes a cut-off with dependence on scales k ^ /co- 
in this paper we will consider a class of models where DM is a mixture of cold and 
warm (with velocity dispersion in the form (1)) components. We define the WDM 
fraction /wdm as 

^ l^wDM n _ n I n I n (a] 

J WDM — Q ; ^^M — 'CDM I ^ 'WDM T i \^ ) 

i 'm 

(note the presence of the baryonic fraction Dr in the denominator). The transfer 
function of ACWDM models is qualitatively different. The characteristic behavior 
of the ACWDM power spectrum as compared to the pure ACDM one for the same 
cosmological parameters is shown on Figs. 1, 2, and can be summarized as follows: 

- The scale where the transfer function departs from 1 depends on the mass m^Rp 
of the WDM component, but not on /wdm? as Figs. 1 and 2 demonstrate. 

- For large /c's, the transfer function approaches a constant plateau. The height of 
this plateau is determined solely by /wdm, as Fig. 2 demonstrates. 



^Eq.(l) is only an approximation to the real sterile neutrino DM spectrum, which can be com- 
puted only numerically (see [(S2] for details). 
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Fi gure 1: Transfer functions Tijt) = [-Pacwdm(^)/-Pacdm(^)] ^ for ''^nrp — 

2 keV 

and different WDM fractions: from left to right, f2wDM/(^wDM + ^cdm) = 
1, 0.8, 0.6, 0.4, 0.2, 0.1. Other parameters are fixed to = 0.05, = 0.3, Vt^ = 0.7, 
h = 0.7. 



2.1 Characteristic scales 

At any given time, the comoving free-streaming wavenumber of the warm component 
is defined as 

fc,.^V^^ = /||| (5) 

where p is the total density of the universe at the time considered, {v) is the velocity 
dispersion of the warm component only (not to be confused with the velocity dis- 
persion of the total cold plus warm component, which depends on /wdm), a is the 
scale factor and H the Hubble rate. This definition is analogous to that of the Jeans 
scale for a fluid. The wavenumber fcpg appears in the perturbation equations, and 
corresponds to the wavelength below which perturbations cannot experience gravi- 
tational collapse due to their velocity dispersion. In the radiation era, when the fluid 
is ultra- relativistic, (v) ~ 1 (using units where c = 1) and fcps evolves like the co- 
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Figure 2: The height of plateau in the transfer function T{k) = 

[PAcwDM(^)/-fAcDM(^)]^^^ depends only on /wdm only and does not change with the 
mass. In these examples, the mass is equal to 2 eV (solid line) or 20 eV (dashed- 
dotted line), while ^wdm/ (^wdm+^^odm) = 0.1 (red) or 0.05 (green). Other parameters 
are fixed to Qb = 0.05, = 0.3, Qa = 0.7, h = 0.7. 

moving Hubble scale, kpg = yJ'ijlaH oc After the non-relativistic transition, 

(f) = {p)/m scales like a~^; with our choice for /(p), the average velocity is given 
approximately by {v) = ?)Ty/m^^p. Hence, during the rest of radiation domination, 
the comoving free-streaming scale remains constant with 

where T° is the current neutrino temperature, and we assumed the scale factor to be 
normalized to unity today. Later on, k^s increases like t^^^ during matter domination, 
and even faster during A domination. Today, assuming the universe to be flat, the 
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free-streaming scale is equal to 



3T0 



1 rj ^NRF 



X 10' 



1 keV/ 



) h/Mpc 



(7) 



Below this scale (for k > /cps), perturbations can never experience gravitational 
clustering. Hence, this scale marks roughly the beginning of the region for which 
the free-streaming effect is maximal and the CWDM power spectrum is reduced by 
a constant amplitude (the plateau in Fig. 2). 

By analogy with the case of light active neutrinos, one could expect that the 
minimum of the free-streaming scale /c™™ (Eq. (6)) would give the point at which the 
ACWDM power spectrum starts to differ from the ACDM one. However, /c™™ does 
not give the right answer, as can be checked on figure Fig. 2. In fact, the largest scale 
affected by free-streaming is nothing but the present value of the particle horizon of 
warm particles with a typical velocity (v), that we will call the (comoving) free- 
streaming horizon ApgH 



A,; 



to 



dt 



am 



da . 



(8) 



This scale can be computed easily if we neglect the impact of A (it is easy to check 
that the contribution of the A-dominated stage to the above integral is completely 
negligible; actually, neglecting also the matter-dominated stage would be sufficient 
for an order of magnitude estimate). We can decompose the integral as: 



Ac,« 



da 3T° r da 

a^H jtinrp Ja^_^^ a 



(9) 



with = T°/mNRP being the scale factor at the time of the non-relativistic transi- 



tion. Using the Friedman equation, this becomes 

1 , 3anr 



Ap 



da 



JO 



da 



J anr 



y/l + tteq/a 



(10) 



where acq = Q^/Qm is the scale factor at the time of matter /radiation equality. In 
the limit Onr <^ 1, the result is simply 



3^0 



FSH 



1 + 6 Arcsinh^ / — ^ 



which corresponds to the wavenumber 



1 + 31n 



^Oeq 




(12) 



10 



Taking ^IrH'^ = 4 x 10 ^ we finally obtain 



"'FSH 



0.1 



flMh^ J Vl kev) 



VMpc . (13) 



One can check in Fig. 2 that this corresponds indeed to the scale at which the 
ACWDM power spectrum starts to differ from the ACDM one. Note that the dis- 
tinction between the free-streaming horizon and the free-streaming scale is not com- 
monly used, because in the case of active neutrinos (which become non-relativistic 
during matter domination) one has k^^-^ ~ in this case, it is not necessary to in- 
troduce the quantity fcpsn, and the scale marking the beginning of the step is usually 
said to be /c™" = k^l- In the case of WDM or CWDM, the distinction is important, 
because between t^r and teq the free-streaming scale Aps remains constant, while the 
free-streaming horizon Apsn increases typically by one order of magnitude. 



2.2 Plateau in CWDM 

We can try to estimate the amplitude Tpiateau of the small-scale power spectrum 
suppression, 

for k » = 8 X 10^ f h/Mpc . (14) 

VI kev / 

Let us consider a ACDM and a ACWDM model sharing the same cosmological 
parameters, and in particular the same Q-om and Q^. For simplicity, we neglect 
the mass of active neutrinos throughout this work. Until the time at which the 
warm component of the ACWDM model becomes non-relativistic, these models are 
indistinguishable. After this transition, but still before matter/radiation equality, the 
evolution of the gravitational potential is fixed by that of the photon and neutrino 
components, which dominate over dark matter; hence, the (slow) growth of cold dark 
matter inhomogeneities, (5cdm = ^Pcdm/Pcdm, is independent of /wdm- We conclude 
that for any wavenumber k, ScDuiO'eq) is the same in the two cases /wdm = and 

/wdm 7^ 0. 

After equality, the evolution of the gravitational potential is fixed by the clustering 
of the dark matter and baryon components. Inside the Hubble radius, the growth of 
5cDM can be inferred from the Newtonian equation of evolution 

r , ^ h H^2f f ^\ (PCDM ~l" Pb)^CDM ~I" PwDm'^WDM r\ 

cicDM + -^CDM = 47rc^a dpTox = o - (15) 

a 2 \ay p^oT 

where we used the fact that 5cdm = soon after equality, and a dot denotes a 
derivative with respect to conformal time. In the ACDM model (with Pwdm = 0), 



T, 



ACWDM 



plateau 



ACDM 



ik). 
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this equation has a well-known solution leading today to 



5cDM(ao) = ( ) 5cDM(aeq) (16) 



where the constant g{ao) < 1 is determined by the value of Q\ = l—Q^, and describes 
the reduction in the perturbation growth rate during the A dominated epoch. It is 
roughly approximated by 

g{ao) - (1 + 0m3{njn^)^/') . (17) 

(see e.g. [98] for details). A more precise numerical evaluation gives g{ao) = 0.779 
for = 0.7. The product ag{a) is usually called the linear growth factor. 

In the ACWDM model, PwDM cannot be neglected in the term Ptot appearing in 
equation (15), but the warm component never clusters on small scales with k > k^g, 
so can be set to zero. Then, a factor 

PCDM + Pb ^^g^ 

Ptot 

appears in the right hand-side of equation (15). During matter domination, this 
factor is equal to (1 — /wdm) and the linear growth factor can be found exactly (like 
in the case of ACHDM models with massive active neutrinos): 

\ l-(3/5)/wDM 

) ScDuiaeq) ■ (19) 

During A domination, there is no simple analytic result. In the case of ACHDM, it 
has been shown in Ref. [98] that a very good approximation is found by just inserting 
the function g[a) inside the parenthesis. In our case, the same ansatz would give 

ScDuiao) = I — — 5cDM(acq) (20) 

V '^cq / 

leading to a reduction in clustering with respect to the ACDM model equal to 

St7='{ao) V «eq ; 
The power spectrum P{k) is defined as the variance of the total matter perturbation 

6p Pb^B + PCDM^CDM + PwDM^WDM 



P Pb + PCDM + Pv 



(22) 
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Today, Sb = ^cdm, and in the ACWDM model (5wdm ^ Scou for k > k^^. So, on such 
small scales, on has 

-Pacdm(^) = (^cdm) 5 -fAcwDM(^) — /wdm) ('^cdm) • (23) 

Using equations (21) and (23), the suppression factor defined in eq. (14) is found to 
be: 

faog{ao) \ ~(^/^)-/^^™" 

^plateau ~ (1 /wdm) I 1 ■ (^4) 

V '-''cq / 

This equation is slightly different from the equivalent result in ACWDM models, 
due to the fact that the WDM component becomes non-relativistic during radiation 
domination. 

Equation (24) encodes the right dependence of Tpiatoau on /wdm, but cannot fit 
accurately the numerical results of a Boltzmann code for two reasons: 

- first, in the above solution, we treated the transition from radiation to matter 
domination as an instantaneous process. In fact, matter comes to dominate in a 
progressive way, and 5cdm ('^) starts to depart from 5cDM'~°(a) slightly before the 
time of equality. Hence, the ratio Sqom [o-eq) / Scdu~^ {o-eq) is actually smaller than 
one and decreases when /wdm increases. This means that equation (24) slightly 
underestimates the effect of WDM. Empirically, we found that this can be corrected 
by increasing the exponential factor from 3/5 to 3/4 (at least, as long as flu 
does not depart too much from 0.3). In figure 3, we compare the result of a full 
Boltzmann code calculation (using CAMB [99]) with the empirical fit 

^plateau ~ (1 ~ /wdm) ( j (25) 

and find very good agreement. In this example, 

for Tomb = 2.726 K, = 0.7, = 0.3, h = 0.7. 

- second, we assumed that in all cases one has 6^ — 5cdm soon after equality. How- 
ever, we are interested here in very small scales {k ~ 100/i/Mpc) for which the 
baryon and CDM overdensities do not converge towards each other before the end 
of matter domination. Hence, the approximations (5cdm oc a or Scdm oc a^-i^/^)f™ 
only apply to a short period: at the beginning of matter domination, baryons do 
not contribute to the right-hand side of equation (15), and the growth of 5cdm 
is slower. Since this effect (governed by the parameter ^b/^m) is similar in the 
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two cases /wdm 7^ and /wdm = 0, one does not expect a strong dependence of 
^plateau on the baryon abundance. As can be seen in figure 3, tfie Boltzmann code 
calculation shows that even radical variations in fie hardly affect the behavior of 
the suppression factor. We conclude that equation (25) is accurate up to a few 
percent for realistic ACWDM models^. 

For small /wdm, one can expand equation (25) to obtain: 
2 6 

-^plateau ~ 1 ^ 2/wDM ~ ^/wDM log Q! ~ 1 — 14/wdm • (27) 

This expression can be compared with Tpj^^^^^^^^ ^ 1 — 8/;^ obtained in the case of active 
neutrinos [100]. 

Note that equation (25) applies today, but the suppression factor at higher red- 
shift can simply be obtain by replacing ao by a = ao/(l + z); so, it is clear that 
the transfer function T{k) is redshift-dependent, with a smaller step-like suppression 
at large z. The function T{k) is only expected to vary with z on scales such that 

<^wdm(-^) < "^cdm 

(z), i.e. for k > kpii{z). This redshift dependence is illustrated in 

figure 4. 



3 hy-a method 

In order to constrain cosmological models, it is always preferable to rely on data 
dealing with linear scales. However, we have seen in the previous section that the 
difference between ACDM and ACWDM (or AWDM) models appears only on small 
scales (typically. A; > 0.1 h/Mpc) which cannot be probed by CMB experiments, or 
by the reconstruction of the galaxy power spectrum in the range in which the light-to- 
mass bias can be treated as scale-independent. Hence, at the moment, the only way to 
discriminate between ACDM and ACWDM (or AWDM) models is to use observations 
of Lyman-a (Ly-a) absorption in the spectrum of distant quasars (QSOs), which can 
be used as a tracer of cosmological fluctuations on scales /c ~ (0.1 — 10) ZiMpc"^, at 
redshifts 2; ~ (2 — 4) where the cosmological perturbations are already affected by 
the non-linear evolution - although not as much as today. 

The Ly-a forest is a dense set of absorption lines observed in QSO spectra. It 
is well established by analytical calculations and hydrodynamical simulations that 
these absorption lines correspond to Ly-« absorption (Is ^ 2p) by clouds of neutral 
hydrogen at different redshifts along the line of sight. This neutral hydrogen is 
part of the warm (~ 10^ K) and photoionized intergalactic medium (IGM). Opacity 
fluctuations in the spectra arise from fluctuations in the neutral hydrogen density, 

"^This statement is meant for the plateau scales k ^ 100/i/Mpc. At much larger k, the situation 
becomes more complicated, since according to the linear theory one has Sb < Sqdm even today. 
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Figure 3: Suppression factor Tpiateau = [-fAcwDM(^)/-PAcDM(^)]"'^ computed today at 
k = 10^ h/Mpc with the Boltzmann code CAMB, for different values of /wdm = 
Qwdm/^m and Qjj. Other parameters are kept fixed: fl\ = 0.7, flu = 0.3, h = 
0.7. The dashed curve corresponds to the analytic prediction of Eq. (25) with a = 
aog{ao)/aeq = 2.75 x 10^. 



from which it is possible to infer fluctuations in the total matter distribution [101, 
102, 103]. For each quasar, the observed spectrum I{z) can be expanded in (one- 
dimensional) Fourier space. The expectation value of the squared Fourier spectrum 
is called the flux power spectrum. The data provide an estimate of the flux power 
spectrum Pp for various redshifts 2 ~ 2 — 4 (for which light rays of wavelength 
Xobs = (1 + 2)1265A pass through the optical window of the Earth's atmosphere). At 
these redshifts, the density perturbations of scales ~ (0.1 — 10) /i Mpc~^ already 
entered into a mildly non-linear stage of gravitational collapse {6p/p > 1). 

The linear matter power spectrum P{k) and the flux power spectrum Ppik) 
are complicated functions of the cosmological parameters. In the rather narrow 
range of scales probed by Ly-a data, the effects of an admixture of WDM might 
be compensated by a change in other cosmological parameters (ag, fluh"^, ng, etc.). 
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Figure 4: Variation of T{k) = [Pacwdm(^)/-Pacdm(^)]^^^ with redshift, for various 
values of m and flwou/ (^wdm + ^cdm); other parameters are fixed to the same values 
as in previous figures. (Left) For m = 0.5 keV, the variation is clearly visible at large 
k > k-ps{z). (Right) For m = 5 keV, the free-streaming wavenumber is too large for 
any variation to be seen on the scales displayed here. 

or even in extra astrophysical parameters in the case of the flux power spectrum. 
Therefore, it is important to perform collective fits to the Ly-a data and other data 
sets (e.g. CMB). Cosmological parameter extraction from combined cosmological 
data sets can be conveniently performed with a Monte-Carlo Markov Chain (MCMC) 
technique, using e.g. the public code CoSMoMC [10 1]). 

In order to predict Pp^k, z) for a given cosmological model, it is necessary to 
perform hydrodynamical simulations (while for CMB experiments probing the linear 
matter power spectrum, it is sufficient to compute the evolution of linear cosmologi- 
cal perturbations using a Boltzmann code like camb). Hydrodynamical simulations 
are necessary both for simulating the non-linear stage of structure formation and 
for computing the evolution of thermodynamical quantities (in order to relate the 
non-linear matter power spectrum to the observable flux power spectrum). The 
characteristic number of samples required for a reliable determination of parameter 
probabilities ranges from 10^ for the simplest cosmological scenario to 10^ (or even 
10^) when more parameters and freedom are introduced in the model. In principle, in 
order to fit Lyman-a data, one should perform a full hydrodynamical simulation for 
each Monte-Carlo sample. This is computationally prohibitive, and various simpli- 
fying approximations have been proposed [105, 50, 61, 60, 65, 106, 62, 107]. Below, 
we discuss the two datasets used in this paper, for which different strategies have 
been implemented. 



16 



3.1 VHS data 



The VHS dataset was compiled in Ref. [95] using 57 QSO spectra from LUQAS 
[108] at z ~ 2.1, and from Ref. [109] at 2; ~ 2.7. Each spectrum was observed 
with high resolution and high signal-to-noise, resulting in clean and robust mea- 
surements. However, this dataset is based on a relatively small number of spectra, 
leading to a large statistical uncertainty. Systematic errors were estimated in a con- 
servative way, and the flux power spectrum was only used in the range 0.003 s/km < 
k < 0.03 s/km (roughly corresponding to scales 0.3 h/Mpc < A; < 3 h/Mpc): at 
larger scales, the errors due to uncertainties in fitting a continuum (i.e. in remov- 
ing the long wavelength dependence of the spectrum emitted by each QSO) become 
very large, while at smaller scales, the contribution of metal absorption systems be- 
comes dominant. It was shown in Ref. [95] that the dependence of the bias function 
b{k, z) = Ppik, z)/P{k, z) on cosmological parameters can be neglected for this data 
set, given the large systematic plus statistical errorbars in VHS, and the limited 
range in /c-space and 2;-space probed by this data [50] . Hence, this bias function can 
be computed for a fiducial model, and the data points Ppiki.Zj) can be translated 
in linear spectrum measurements P{ki,Zj). This approach is implemented in the 
pubhc module lya.f90 [65] distributed with CoSMoMC. We will use this dataset 
in section 5 to derive conservative estimates. 

3.2 SDSS data 

For the Ly-a data obtained by the SDSS collaboration [96], a different approach 
has been implemented in Ref. [61]. The data points were obtained from 3035 quasar 
spectra with low resolution and low signal-to-noise, spanning a wide range of redshifts 
{z = 2.2 — 4.2). Dealing with low resolution and low signal-to-noise QSO spectra and 
extracting the flux power is a very difficult task as shown in [61]. We refer to their 
analysis for a comprehensive study of continuum fluctuations removal, metal lines 
contamination, presence of damped Ly-a systems, resolution of the spectrograph 
and noise level in each redshift bins. All these effects need to be properly taken into 
account, and not modelling them properly would impact the obtained flux power in a 
way that is difficult to predict. In the following, we will use the flux power provided 
by the SDSS collaboration, introducing nuisance parameters for the resolution and 
noise in each redshift bin as suggest by [96, 61], and implicitly assuming that all the 
contaminants above have been either removed or properly modelled in their analysis. 

The low resolution of SDSS data implies that the flux power spectrum cannot 
be measured on small scales (with k > 0.02 s/km). However, on larger scales, the 
low resolution and signal-to-noise are compensated by the statistics: with such a 
large number of quasars, statistical errorbars on the flux power spectrum are much 
smaller than for VHS data. Hence, using a fixed bias function would be inappropri- 
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ate. Instead, Ref. [96, 61] used a large number of fast hydro-particle- mesh (HPM) 
simulations [110, 106] densely sampling the cosmological parameter space, calibrated 
from a small number of full hydrodynamical simulations. These HPM runs could 
produce an inaccurate flux spectrum estimation for cosmological models departing 
from the fiducial points where the hydrodynamical simulations were performed, and 
a more reliable calibration is probably needed. 

A different approximation scheme was used in Ref. [62] using the same data. 
For each wavenumber and redshift, the flux power Pp{k,z) was Taylor-expanded 
around a fiducial model at first order in cosmological parameters, using of the order 
of full hydrodynamical simulations (here ~ 20 is the number of relevant cos- 
mological/astrophysical parameters). This method can also become inaccurate far 
from the fiducial model, especially when the likelihood and/or parameter posterior 
probability deviates significantly from a multivariate Gaussian distribution. How- 
ever, in both methods, it can be checked a posteriori that the inaccuracy introduced 
by interpolating/extrapolating around the points where hydrodynamical simulations 
were performed is negligible with respect to data error bars, at least in the range of 
cosmological models producing reasonable fits to the data. 

We deal with uncertainties related to Ly-a physics in the same way as in Ref. [61, 
62]. Namely, we use 9 "astrophysical" nuisance parameters. Two describe the ef- 
fective evolution of the optical depth, six describe the evolution of the parameters 
7 and Tq entering the relation between the temperature of the intergalactic medium 
and its overdensity, T = Tq{1 + 5)'^"^. The parameter Tq was allowed to vary in the 
range < Tq < 10^ K, while < 7 < 2. Note that the range of parameters chosen is 
very wide and embraces also the recent results of [55] that found evidences of an in- 
verted (7 < 1) temperature-density relation. The parameters describing the thermal 
evolution are modelled as a broken power-law at 2 = 3 as ?/ = A{^-^)^, with one 
parameter for the amplitude and two for the slopes at 2; < 3 and 2; > 3 (the slopes 
are left free to vary in the range [—6, 6]). This is done in order to better mimic more 
complex evolutions able for example to reproduce the Helium II reionization at z ~ 3 
as suggested by SDSS observations [111]. The evolution of the effective optical depth 
is modeled as a single power-law (two parameters: one for the amplitude and one 
for the slope) within a very wide observational range. Another parameter describes 
the contribution of strong absorption systems (Damped Lyman-a systems) and this 
is left free following exactly the suggestions of Ref. [61] and applying the correction 
of Ref. [112]. 

Finally, we have a total of 9 additional parameters that model uncertainties in 
the correction to the data - eight parameters for the noise correction in various red- 
shift bins, and one parameter describing the error of the resolution correction (these 
parameters do have strong a-posteriori priors as in Ref. [61]). All these parameters 
as well as those describing corrections for damped/high column density systems are 
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constrained as suggested in [112] and described in [62]. 

When deriving bounds on the parameters of interest below, we marginalized over 
these 18 parameters, and the final estimates of the astrophysical parameters are 
within the observational bounds. 

In order to take into account the contribution of Silll metal lines to the power 
spectrum, we followed an approach similar to the one suggested in [61]. Namely, we 
introduced an extra parameter, describing the contribution of Silll lines to the power 
spectrum, and fixed it to the maximum correction suggested by the SDSS collabo- 
ration [61], motivated by the following arguments: i) the metal lines contribution 
peaks at smaller scales than those investigated here; ii) at large scales the effect is 
an overall suppression of power which is already taken into account by the fact that 
we are considering a wide range of possible mean effective optical depths ([108]); iii) 
the effect of modelling further the Silll metal contribution by allowing ^-dependence 
has been found to have a small impact on the final results by Ref. [61]. 

Within the range allowed by the data, the effect of ACDM cosmological parame- 
ters and astrophysical parameters on the flux power spectrum is nearly linear, as can 
be checked a posteriori from the fact that all marginalized probabilities are almost 
Gaussian. This justifies the use of a first-order Taylor expansion. The authors of 
Ref. [62] checked this explicitly for some parameters, by performing a second order 
Taylor expansion or by trying to estimate the flux power in regions of the parameter 
space that are very far from the best-fit with a full hydro dynamical simulation. They 
usually found good agreement between the various methods at the percent level in 
the simulated flux power. 

In the present work, we want to obtain constraints on mixed ACWDM mod- 
els. The parameter probability distribution is strongly non-Gaussian with respect to 
the extra parameters characterizing this set-up: the density fraction of warm dark 
matter, and its mass (or velocity dispersion). Hence, a simple first-order Taylor 
expansion would lead to inaccurate results. In section 6, we extend the method of 
[62] by interpolating between a grid of fiducial models in this two-dimensional sub- 
space, while keeping a simple Taylor expansion in other parameters. To this end, we 
performed a grid of 8 additional simulations (for two values of the masses and four 
values of the warm dark matter fraction). In the next section, we will discuss the 
robustness of these simulations, especially as far as the issue of initial velocities is 
concerned. 

4 Thermal velocities in ICs 

In N-body or hydrodynamical simulations, the difference between Initial Conditions 
(ICs) for CDM and WDM does not reside entirely in a modification of the power 
spectrum, but also in non- negligible thermal velocities for the warm component. For 
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Thermal Relics (TR) assumed to decouple from thermal equilibrium when they are 
still relativistic, and characterized by a Fermi-Dirac distribution with temperature 
Ttr, the average velocity is given by 

^ 3.151 r„(.) _ + n_k^\ /T^^ 

^ ' m„ V 100 y V '"t» J \1KJ ^ ' 

However, in this paper, we assume that WDM consists in sterile neutrinos produced 
through a Non-Resonant Production (NRP) mechanism, namely, non-resonant oscil- 
lations with active neutrinos [15, 84, 82]. In this case the distribution can be roughly 
approximated by the renormalized Fermi-Dirac distribution of Eq. (1), corresponding 
to an average velocity 

^ 3.151T,(z) ^ 3.151(4/11)^3(1 + ;,)TcMB ^ /l±f ^ /^ l keV 

mNRP mNRP V 100 / V "^NRP 

(29) 

If the mass and the relic density I^wdm^^ are specified, the parameter Ttr (for TR) 
or X (for NRP) can be inferred from (Ttr/T,^)^ = X = ^wdm^^(94 eV/m). Thermal 
relics and sterile neutrinos sharing the same mass and density correspond to different 
initial velocities. For instance, assuming that Q^dm^^ = 0.12, we find the following 
velocity dispersion at z = 99 in each of the two cases: 

m = 10 keV (^nrp) ~ 1-6 km/sec (i'tr) ~ 0.16 km/sec 

m = 5 keV (^nrp) ~ 3.1 km/sec {vtr) ~ 0.41 km/sec (30) 

m = 0.5 keV (wnrp) ~ 31 km/sec ("Wtr) ~ 8.8 km/sec 

These velocities should be compared with the peculiar velocities which particles 
acquire when clustering. In most codes, the latter are initialized using the so-called 
Zel'dovich prescription [113, 114, 115]: so, later, we will make the distinction between 
thermal velocities and Zel'dovich velocities. Whenever the thermal velocities are 
negligible with respect to the Zel'dovich ones, the final result is expected to be 
insensitive to the proper implementation of thermal velocities, which can then be 
omitted from ICs. In a typical simulation, the Zel'dovich velocities are of below or 
around 20 km/s at z = 99. Hence, for the NRP case, we roughly expect. Hence, 
for the NRP case, we roughly expect that thermal velocities are unimportant for 
m > 5 keV. 

This can be checked explicitly. Our numerical simulations are performed with 
the GADGET-ll code [116], extended in order to compute the Ly-a flux power spec- 
trum [60]. First, we address the case of pure AWDM models with Qcdm = 0. We fix 
initial conditions a.t z = 99. The initial power spectrum is computed with camb, 
and for each warm particle a random thermal velocity is eventually added to the 



j 15.7km.s"^ 



20 



Zel'dovich velocity. Thermal velocities are distributed with a Fermi-Dirac probabil- 
ity. 

Notice, however, that a particle in the simulated ICs represents a region with size 
/ = Lhox/N^^^ (where is the number of particles) . Therefore, a huge amount of DM 
particles are represented as one body in N-body simulations. Formally, the velocity of 
such a collection of particles is zero, as the velocities of real DM particles are random 
and they average to zero. The assignment of a Fermi-Dirac thermal velocity to such a 
simulation particle is justified provided that / is smaller than the smallest scale probed 
by our data set, and larger than the extra distance traveled by the particle, due to 
thermal velocities, when it moves along Zel'dovich trajectories. In our simulation, 
we take = 400^ = 6.4 x 10^ particles in a box of comoving size L^ox = 60Mpc/h, so 
each resolution element in the ICs measures / ~ O.lSMpc in comoving units, while the 
smallest scale probed by our Ly-a flux power spectrum corresponds to a comoving 
wavelength of the order of 1 Mpc. In addition, the extra distance travelled by the 
particle can be estimated conservatively by multiplying the velocity f ~ 10 km/s by 
the amount of time between 2; = 99 and 2 = 2, t ~ 3 Gyr: this gives 0.04 Mpc,which 
is smaller than the resolution element probed by the particle in the ICs. Hence, our 
approach is justified. A similar way to include velocities in the WDM simulations 
was used e.g. in [117]. Note that the simulated flux power has been corrected for 
resolution using exactly the same procedure described in [68], so even if most of the 
simulations have been ran at somewhat low resolution we corrected for it using the 
higher resolution runs. 

In Fig. 5, we show the ratio of the AWDM flux power spectrum to the ACDM 
one at k = 0.14 h/Mpc and z = 2.2 for 3 masses (0.5,5, lOkeV), with and without 
thermal velocities in ICs. One can see that for masses above ~ 5 keV, the influence 
of velocities is below the 1% level, which is smaller than the error bars of the SDSS 
data points (in case of the VHS dataset, the influence of velocities is unimportant 
also for smaller masses due to the larger error bars). For smaller masses, including 
thermal velocities in the initial conditions becomes essential. 

This suggest that also in the case of a mixed ACWDM, we can perform simulations 
without thermal velocities provided that the mass of the WDM NRP component is 
above ~ 5 keV. This simplifies considerably the issue of initial conditions, because 
in this case we can treat the mixed CDM-I-WDM components as a single cold fluid 
with negligible thermal velocities, and a proper initial power spectrum (equal to the 
total linear matter power spectrum at high redshift, as computed by camb). 

5 Results for VHS data 

We performed a Bayesian parameter estimation for the ACWDM model, using the 
following data set: VHS Ly-a data; WMAP5 [118] and small-scale CMB experiments 
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Figure 5: Dependence of the flux PS of AWDM models (divided by that of ACDM) on 
thermal velocities in hydrodynamical simulations. The points show p^WDM^'pACDM 
at /c = 0.14 /i/Mpc and z = 2.2, extracted from GADGET-II simulations for three 
different masses, with and without thermal velocities in the initial conditions. The 
lines are the results of linear interpolation between these points. 



(ACBAR [119], CBI [120], Boomerang [121, 122, 123, 124]); galaxy power spectrum 
from the SDSS LRGs [125]; supernovae from the SNLS [126]. We take fiat priors on 
six cosmological parameters describing the ACDM model {Qumh'^, ^^b^^, t, Us, As)^, 
on two extra parameters describing the warm sector (see below), and we marginalize 
over nuisance parameters associated with the data sets, in the way implemented in 
the public version of CoSMoMC (compiled with the file lya.fQO). The two extra 
parameters are chosen to be: 



first, the warm fraction of dark matter. 



^CDM ~l~ 



(31) 



^See e.g. [104] and http://cosmologist.info/cosmomc/readme. html for explanation of this 
parameter choice. 
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Figure 6: Preferred region for pairs of cosmological parameters in the ACWDM 
model, using VHS, WMAP5 and other cosmological data sets. The black lines repre- 
sent the 68%CL and 95%CL contours of the two-dimensional probability marginal- 
ized over other parameters (as defined in the Bayesian approach), while the colors 
show the likelihood averaged over other parameters. The first five parameters Q^h"^, 
(w)wDM, ns have flat priors over the displayed range. Our analysis also 
assumes fiat priors on the logarithm of the primordial spectrum amplitude and on 
the angular diameter of the sound horizon, but instead of these parameters we show 
here erg and Hq, which have a more straightforward interpretation. In this plot, we 
do not show the reionization optical depth r and all the nuisance parameters enter- 
ing into the analysis. The one- dimensional marginalized probabilities are displayed 
along the diagonal. 23 
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Figure 7: (Left) Preferred region for the parameters (F^dm, {'v)wdm) of the ACWDM 
model, using VHS, WMAP5 and other cosmological data sets. The black lines repre- 
sent the 68%CL and 95%CL contours of the two-dimensional probability marginal- 
ized over other parameters (as defined in the Bayesian approach), while the colors 
show the likelihood averaged over other parameters. (Right) Same plot with the hor- 
izontal axis translated in terms of {IkeV/m^^p) = (f)wDM/(157 m/s), see equation 
(29). 



which differs from the quantity /wdm defined in Eq. (4) by (^om/^a/) since 
the denominator accounts for dark matter instead of total matter; indeed, 
in the analytical solutions of section 2, /wdm is the quantity which appears 
naturally; in a parameter extraction, -F^dm is more convenient because of its 
clear interpretation: it varies from Fwdm = for pure cold to Fwdm = 1 for 
pure warm dark matter. Instead, /wdm = 1 would describe a model with 
negligible CDM and baryonic components, providing very bad fits to the data, 
and an eventual upper limit on /wdm would not be straightforward to interpret 
physically. 

• the average velocity of warm particles today in km/s: (t')wDM- This parameter 
is convenient because, unlike the particle mass, it has a direct physical effect on 
the power spectrum; hence our results apply to both NRP and TR models (and 
presumably to a much wider class of WDM models). Bounds on (u)wdm can be 
translated immediately in terms of rrinRp using Eqs. (29); limits on are not 
so straightforward to obtain since they depends on the posterior distribution 

of ('y)wDM) -^WDM and l^DM- 

Note that varying i^wDM and (f)wDM is equivalent to tuning respectively the size and 
the location of the step-like suppression in the matter power spectrum (at least at 
the linear level), as shown in Section 2. 
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The main goal of the various MCMC runs performed in this work is to derive 
some (Bayesian) credible intervals*^ for the parameter (u)wdm, starting from a flat 
prior on this parameter in the range < (w)wdm < oc. These results can be imme- 
diately translated in credible intervals for m'^p, using the relation (1 keV)/mNRp = 
('i^)wDM/(0.157 km/s), see equation (29) with z = 0. It would not be convenient to 
derive directly some credible intervals on mNRP itself, since the ACDM limit (which 
is perfectly allowed by the data) corresponds to However, each 95% 

Confidence Level (CL) upper bound on (f )wdm can be expressed as a lower bound on 
the mass, that we will denote as m^^p. We warn the reader that this bound should be 
interpreted with care, since we are not using a fiat prior on the mass, but rather on 
the velocity, i.e. on the inverse mass. Hence the quantity m^^p represents a Bayesian 
95% CL upper limit on the velocity expressed in a useful way. 

We first ran COSMOMC with -Fwdm fixed to one, in order to update bounds on 
the pure AWDM model. We found a credible interval < (t')wDM < 0.10 km/s 
(95% CL), which can be translated into the lower limit m^^^p = 1.6 keV. In refer- 
ence [65] (which also used VHS Ly-a data, but WMAPl instead of WMAP5), a 
stronger bound m^^p = 2.0 keV was found. This is surprising at first sight, since 
the CMB temperature and polarization spectra are better measured by WMAP5. In 
addition, WMAP5 favors the lowest range of rig and ag values allowed by WMAPl, 
so one could expect that a cut-off in the power spectrum due to WDM is even more 
unlikely, leading to stronger bounds. Actually, the constraints on this cut-off come 
entirely from the VHS data: so, for AWDM models, the real expectation is that the 
mass bounds should remain almost insensitive to variations in the CMB data. Still, 
a possible explanation of this 20% discrepancy is that for AWDM (and ACWDM) 
models, the convergence of MCMCs requires a huge amount of samples. In this 
work, we deliberately accumulated a very large number of chain points, and checked 
thoroughly that our chains are fully converged'. The older result might have been 
affected by poorer chain convergence. 

^See the Appendix A for a comparison of Bayesian credible intervals and frequentist confidence 
limits. 

^ This particular AWDM run is based on 6 independent chains. After accumulating a total of 
^ 30 000 chain points, we reached a convergence criterion (i? — 1) ~ 0.02 (based on the second half 
of the full chains, see Ref. [I<I4] for a definition of i?), and a bound m^lp ~ 2.0 keV. However, it is 
only after accumulating ^ 100 000 chain points that the bound stabilized to TOnrp = 1.6 keV, with 
(i? — 1) ^ 0.007. We know that this bound is accurate because we monitored the evolution of the 
mass bound versus execution time, and checked that its variations remain negligible (less than 5%) 
during the second half of the run; also, we checked that if we vary the fraction of points to remove 
at the beginning of each chain in the range 20-80%, or if we remove one of the chains, the bound 
remains precisely equal to 1.6keV. In the next section on SDSS results, we will use CoSMOMC with 
a more convenient convergence criterion than R in order to establish the precision and stability of 
aU 95% CL credible limits 
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Notice that our 95% CL credible limit on rrij^^p in AWDM models is very close to 
that obtained from phase-space density arguments [21] applied to newly discovered 
dwarf spheroidal satellites of the Milky Way [25], which give mNR? > 1.77 keV (see 
also [34]). 

We then ran CoSMoMC for the full ACWDM model, with one more parameter 
-^wDM- Our results are summarized in Figure 6. The plots on the diagonal show the 
probability of each parameter. The plots off-diagonal show the (Bayesian) credible 
limits (68% CL or 95% CL) on pairs of parameters. 

Figure 7 focuses on the two-dimensional probability distribution in the plane 
(FwDM, {v)wdm) or (FwDM, (lkeV)/mNRp). As F^dm approaches one (pure WDM limit), 
the 95% CL contour reaches (t')wDM ~ O.lOkm/s, which corresponds to a lower bound 
on the mass rri^^p = 1.6 keV (coincidentally, this is identical to the bound derived 
from the one-dimensional probability of (f )wdm in the AWDM case). In the opposite 
limit, as (v)wdm grows, we qualitatively expect that the contour will tend towards 
a constant value, since below a particular value of -Fwdm the impact of the WDM 
component on the linear power spectrum will become too small for being detectable 
by the VHS data, given the error bars and calibration uncertainty of this data. We 
find indeed that the 95%CL contour reaches asymptotically -Fwdm — 0.1: that is, an 
admixture of about 10% of WDM with an arbitrary value of the mass within our 
prior bound {m^^p > 160 eV) is compatible with the VHS Ly-a data*^. 

6 SDSS results 

We now consider the bounds that can be obtained by combining two data sets: 
WMAP5 [118] and the SDSS Ly-a data introduced in section 3.2. We did not 
include any other experimental results in this section. 

6.1 Pure WDM model 

Before applying the method described in section 3.2 to the ACWDM model, we up- 
dated the analysis of [67] for a pure AWDM model. The new ingredients in this work 
are (i) the use of WMAP5 rather than WMAP3 data; (ii) the suggestion of Ref. [55] 
to marginalize the astrophysical parameter 7 over a more conservative range (recall 7 
describes the relation between the local temperature of the intergalactic medium and 
its local overdensity, T = Tq{1 + 6)'^~^). Apart from these new inputs, we performed 
several tests concerning the validity and robustness of the method, trying different 

^Notice, however, that a ACWDM model with 10% of WDM and muRp ~ 160 eV is ruled out by 

phase-space density considerations [25]. Indeed the bound ttinrp > 1.77 keV presented in [25] scales 
1/3 

as /wDM- Therefore when reducing the fraction of WDM from 100% to 10%, it diminishes by 50% 
only and becomes m^-np > 0.89 keV. 
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Figure 8: Dependence of the interpolated flux PS of AWDM models (divided by 
that of ACDM) on thermal velocities in ICs, here for k = 0.018s/ km and z = 4.2. 
The points are obtained with hydrodynamical simulations; the curves correspond to 
a global quadratic flt of these points. 

schemes of interpolation between the points where the hydrodynamical simulations 
were performed, implementing or not thermal velocities in the initial conditions, test- 
ing extensively the convergence of MCMCs, and comparing Bayesian and frequentist 
limits. For these tests we performed GADGET-II simulations for three different 
values of the mass, rriNRP = 0.5, 5, 10 keV (instead of a single one, used in [67]). 

6.1.1 Bayesian credible interval 

It is important to check whether including or not thermal velocities into the ICs 
has an impact on the WDM mass bounds. The analysis of section 4 shows that 
thermal velocities are expected to change signiflcantly the simulated flux PS at least 
for mNRp < 0.5 keV, while the effect of velocities remains very small at least for 
'^NRP ^ 5 keV. Since previous works [66, 67] obtained a (Bayesian) lower bound of 
the order of m^^p = lOkeV (95% CL), one expects the influence of thermal velocities 
in ICs onto the lower mass bound to be small. To analyze this issue in details we 
performed a number of tests. 

In our analysis, the dependence of Pp^°^{k, z) on ttInrp for each value of {k,z) 
is obtained by interpolation, following a global quadratic flt to the four points for 
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Figure 9: Probability distribution for the parameter IkeV /m^^p (which is identical 
to (t')wDM/(157 m/s)) in the pure AWDM case, using WMAP5 and SDSS Ly-a data, 
and starting from a flat prior on this parameter; the solid line shows the marginalized 
posterior, while the dashed line shows for comparison the average likelihood. The 
normalization of the y-axis (the probability) is arbitrary. 

which we have exact results (ttinrp = 10, 5, 0.5 keV, and the pure ACDM 
oo). The result of such a fit is shown in figure 8 for particular values of k and z. 
For mNRP < 5 keV, the flux PS is strongly affected by the presence or absence of 
thermal velocities in the hydrodynamical simulations with m^^p = 0.5 keV, even with 
some level of arbitrariness since various interpolation schemes could be used between 
mNRP = 0.5 keV and mNRP = 5 keV. Even the values of Pp for mNRP > 5 keV are 
slightly affected at the per cent level by the difference between the two interpolations^. 
Thus, the dependence on velocities in ICs can potentially affect the pure WDM result, 
which depends on the flux spectrum variation for mNRP > 10 keV. 

We have performed MCMC runs for both types of simulated flux PS (with and 
without thermal velocities in ICs), using a Bayesian top-hat priors on the parameter 
< {v)wDM < 0.3 km/s (corresponding to mNRp > 0.5 keV). At the 95% CL, we 
found < (f)wDM < 11-6 m/s (with velocities) or < (f)wDM < H-l m/s (without 

^ We have also compared difTerent ways of interpolating the flux PS between the four points at 
our disposal for each (k, z) and found that above jtiwdm ^ 5 keV, the difference between the various 
interpolation schemes is at most of the order of several per cents for each P^™'°^'(fc, z), which is 
comparable with the error bars on the data points of the SDSS set. The resulting errors in mjJIp 
are less than 10%. 
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velocities), corresponding to the bounds m^^p = 13.5 keV (with velocities) or m^^p = 
14.1 keV (without velocities). In figure 9, we show the marginalized probability 
distribution for the variable 1 keV/mNRp, which is identical to (t')wDM/(157 m/s) (see 
equation (29) at 2 = 0). 

As any credible interval limits obtained with an MCMC method, the above values 
of m^^p have some uncertainty related to the fact that there is no clear and unique 
criterion of convergence for MCMCs [104]. It is always difficult to decide when 
the bounds are considered to be stable and when to stop CoSMoMC. To quantify 
this ambiguity, we used some convergence criteria implemented into CoSMoMC: 
MPI_Limit_Converge_Err = 0.3 and MPI_Limit_Converge = 0.025, which mean that 
the run stopped when the 95% CL credible limits were estimated to be accurate 
at the level of 30% of the standard deviation for each parameter. In the case of 
{v)wDM, the standard deviation was found to be cr = 3.6 m/s, so the error on the 
bounds coming form the finite degree of convergence of the chains was of the order 
of ~ 1 m/s. This means that our bounds for the velocity and for the mass are only 
accurate up to ± 10%. Indeed, this refiects the typical variations that one observes 
when varying the fraction of the chains included in the analysis, when removing 
some of the chains, etc. Including this source of error, the two bounds obtained with 
and without thermal velocities are found to be compatible with each other, and are 
summarized by t^nrp — 12.1 keV. This final result for the AWDM case, which takes 
into account the systematic uncertainty of the MCMC method, depends neither on 
the interpolation scheme, nor on the inclusion thermal velocities in ICs. 

This discussion suggests that in the mixed ACWDM case, we should use a top- 
hat prior < (f)wDM < 0.03 km/s (corresponding to ttIwdm > 5 keV) in order to 
obtain conservative bounds; in this way, our final results will not depend on the 
treatment of thermal velocities; and more importantly, they will not be affected 
by the interpolation scheme, which would remain somewhat arbitrary even if the 
thermal velocities were correctly implemented in gadget. In order to obtain robust 
results below mwDM ~ 5 keV, one should perform more hydrodynamical simulations 
for intermediate values of the mass, and perform extensive tests of the convergence 
of GADGET results w.r.t thermal velocity effects. This analysis is beyond the scope 
of the present work. 

6.1.2 Comparison with frequentist approach 

In the previous section, we found that the 95% CL credible interval for (w)wdm cor- 
responds to a lower limit on the mass m^^^ = 13.5 keV (not including the additional 
uncertainty related to MCMC convergence issues, which led us to lower the final 
bound to TTi^^p = 12.1 keV). It is well known that the Bayesian approach finds the 
most probable parameter values given the data (see Appendix A for a short introduc- 
tion to Bayesian and frequentist analyses). Therefore the Bayesian 95% CL limits 
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Table 1: Dependence of the best likelihood in our MCMCs (for SDSS Lyman-a and 
WMAP5 data) on the mass m^^p. The five lines correspond to distinct sets of chains, 
obtained by running CoSMoMC with fixed values of mNRp: CDM limit (equivalent 
to 

^^NRP — ^^^^^ 

), and mNRp equal to Tw^^-p divided by 1, 1.2, 1.5 and 2.0, where Tn^^-p 
is the mass value corresponding to the Bayesian 95% CL upper bound on (w)wdm- 
Notice that these five cases have the same number of degrees of freedom (data points 
minus parameters). 



on each parameter border the region containing 95% of the parameter probability 
(marginalized over other parameters). We are also interested in the answer to a 
different question: whether for a particular value of m^^p there exists at least one 
model which fits the data well enough compared to pure CDM (regardless of how 
many models fit the data equally well). Addressing this question is particularly use- 
ful for excluding conservatively some values of the mass, since it also tells beyond 
which value there is not a single model providing a good fit to the data. 

If the likelihood of the AWDM model was a perfect multivariate Gaussian, we 
would expect that m^|p is the value of the mass for which the maximum likelihood 
(max[£|^^„95^]) is equal to the absolute maximum (max[£]) times e~^. In the general 
case, this is not necessarily the case, and the frequentist bounds can differ significantly 
from the Bayesian ones. 

To estimate the frequentist confidence limits on mNRp, we performed a number of 
CoSMoMC runs with fixed values of the mass. In Table 1, we quote the maximum 
likelihood in each set of chains. Variations in the logarithm —2 log(£bost) (which 
would be x^-distributed in the limit of a multivariate Gaussian) can be interpreted 
as variations in the goodness-of-fit for each value of the mass. By looking at Table 1, 
one sees that the best-fit model for m-NRp = 13.5 keV fits the data almost as well as 
the pure CDM model (the differences Ax^ = 1.2 is not statistically significant). 

As mentioned in Appendix A, the value of the mass for which —2 log(£bcst) wors- 
ens by = 4 (resp. Ax^ = 9) is an approximation for the 95%CL (resp. 99.7%CL) 
lower bound for m^^p, conventionally quoted as the 2a (resp. 3o" bound) in many 
data analyses based on the frequentist approach. In our case, this happens for 
ttInrp — 9.5 keV {2a bound) and m^RP — 8keV (3(j bound) using linear interpolation. 
Hence, mNRP — 8 keV appears to be a robust ?>a lower bound on the WDM mass 
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Figure 10: Dependence of Ax^ = — 21og£ on the WDM mass m^j^p. The values are 
taken from Table 2, and the dashed line linearly interpolates between the points. 

(NRP case) in the pure AWDM case. 

The above derivation of frequentist bounds could be improved by studying more 
values of the mass, finding the maximum likelihood with higher accuracy, and com- 
puting the exact values of Ax^ corresponding to 95% and 99.7% bounds (for non- 
Gaussian distributions, these are not exactly equal to 4 and 9). However, we believe 
that our treatment is sufficient for concluding that values mNRP > 8 keV are allowed 
by present data at the 3a level. 

6.2 SDSS results for CWDM 

In order to analyze the full ACWDM model, we performed extra GADGET-II sim- 
ulations for a grid of points with m^^p = 5, 10 keV and -FWdm = 0.2, 0.5, 0.9. We 
have seen that for such values of the mass, thermal velocities in the ICs can be 
neglected. Hence, we treated the cold plus hot mixture as a single species with 
no thermal velocities, distributed initially according to the appropriate linear power 
spectrum computed by CAMB at 2; = 99. Including ACDM and AWDM models on 
the edges, we have a total grid of 3 x 5 models with mNRP = 5 keV, 10 keV, 00 and 
FwDM = 0,0.2,0.5,0.9,1. For each pair {k,z), we find the value of P^cwdm ^^^^ 
given point (mNRp, -^wdm) by bilinearly interpolating within this grid. 
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Figure 11: Preferred region for pairs of cosmological parameters in the ACWDM 
model, using SDSS Lj-a and WMAP5 data sets. The black lines represent the 
68%CL and 95%CL contours of the two-dimensional probability marginalized over 
other parameters (as defined in the Bayesian approach), while the colors show the 
likelihood averaged over other parameters. The five parameters l^e/i^, ^dm^^, ng, 
FwDM, {v)wDM have fiat priors over the displayed range. Our analysis also assumes 
fiat priors on the logarithm of the primordial spectrum amplitude and on the angu- 
lar diameter of the sound horizon, but instead of these parameters we show here Ug 
and Hq, which have a more straightforward interpretation. In this plot, we do not 
show the reionization optical depth r and all the astrophysical or nuisance param- 
eters entering into the analysis. The one^^mensional marginalized probabilities are 
displayed along the diagonal. 




Figure 12: (Left) Preferred region for the parameters {Fwdm, {v)wdm 

) of the ACWDM 

model, using SDSS Lj-a and WMAP5 data sets. The black lines represent the 
68%CL and 95%CL contours of the two-dimensional probability marginalized over 
other parameters (as defined in the Bayesian approach), while the colors show the 
likelihood averaged over other parameters. (Right) Same plot with the horizontal 
axis translated in terms of (IkeV/mNRp) = (t')wDM/(157 m/s), see equation (29). 



6.2.1 Bayesian credible region 

Having performed extensive tests for the case of pure AWDM models, we applied the 
same type of analysis to the full ACWDM model. We ran CoSMoMC varying all 
parameters, including < (f)wDM < 0.03 km/s (corresponding to mj^^p > 5keV) and 
< -^wDM ^ 1- The resulting two-dimensional joint probabilities are shown in Fig. 11 
and Fig. 12. The data clearly prefers the pure ACDM limit, which corresponds to the 
two axes -Fwdm = and (f )wdm = 0, but significant deviations from this model are still 
perfectly allowed by the data. One can see that as -Pwdm approaches one, the 95%CL 
limit on (F^dm, (1 keV/mNRp)) tends towards {IkeV / m^J^p) = 0.1, i.e. rrij^^p = lOkeV. 
This bound is weaker than the 12.1 keV limit obtained in the pure AWDM analysis, 
because it is derived from a joint two-dimensional probability.^'^ There is a number 
of models with -F^dm < 1 and masses of the order of 5 keV to 10 keV which lie in 
the region preferred by the data. In particular, for F^om lower than 35%, the data is 
compatible with the smallest mass included in our prior range [m^^p = 5 keV). We 
did not include smaller masses in the analysis for the reasons explained at the end 
of section 6.1.1, but it is interesting to see that for mNRP ~ (5 — 7) keV the bound 
on -F^DM reaches a flat asymptote; this suggests that for such a small WDM fraction, 
the amplitude of the suppression in the small-scale flux power spectrum is so small 
that the Ly-a data is compatible with any arbitrary value of the free-streaming scale. 

^"The 95%CL limit in figure 12 contains 95% of the posterior probability in the full 
(^wDM, (1 keV/mNRp)) plane. It is clearly different to look for the restriction of this limit to i^wDM = 1, 
and to compute the range containing 95% of the posterior probability along the i^wDM = 1 axis. 
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Table 2: Values of — log(£best) (which play the role of for runs with fixed Tfij^^ip 

and -FwDM- All runs have the same number of parameters, equal to that of CDM. 

In this case, considerably smaller masses are likely to be allowed, and we would in 
fact need to add data on larger scale (like galaxy redshift surveys) in order to better 
constrain the free-streaming scale and the mass. This would lead to mass limits 
corresponding to hot dark matter rather than warm dark matter particles. 

6.2.2 Comparison with frequentist approach 

Similarly to section 6.1.2 we have performed a number of CoSMoMC runs on a 
grid of points with mNRP = 5, 10 keV and Fwdm = 0.2, 0.5, 0.9, 1.0, varying the same 
parameters as in the pure ACDM case. For each of these runs, we report in Table 2 
the value of — log(£best) in our chains, which plays the role of The dependence 

of on 

FwDM for ''^NRP — 5 keV is presented in Fig. 13. We see that if the mass is 
assumed to take this particular value, the 2a confidence limit (Ax^ = 4) becomes 
FwDM ^ 40%, while the 2a limit (A^^ = 9) allows for about 60% of WDM. We 
can try to compare this case with the results of our Bayesian analysis, for which 
the joint 95% CL limit on (Fwdm, (1 keV/mNRp)) reaches Fwdm = 0.35 in the limit 
"^NRP = 5 keV. In the frequentist approach, the joint 95% CL on a pair of parameters 
is approximated by Ax^ = 6.2. This corresponds to an upper bound -F^dm — 0.50, 
less restrictive than the comparable Bayesian limit. For rrif^^p = lOkeV, all values of 
the WDM fraction (up to pure AWDM) are allowed, in agreement with the results 
of section 6.1.2. 

7 Systematic uncertainties 

The bounds derived previously assume that all systematic uncertainties are properly 
taken into account in the WMAP5, VHS and SDSS Ly-a data likelihoods imple- 
mented in CoSMoMC. Systematic errors are of course hard to estimate. Here, we 
want to summarize the assumptions concerning systematic errors on which our results 
are based, and to describe a number of related tests. 
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Figure 13: Dependence of — log(£best) 

on Fmvdmi assuming that the mass is fixed to 
'^NRP = 5 keV. Horizontal dashed fines sfiow tfie levels Ax^ = 1,4,9, corresponding 
to tfie 1, 2, 3cr confidence limits. 
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The systematic uncertainties entering into our final mass bounds can be roughly 
subdivided into: 

1. Observational uncertainties related to the way the Lj-a forest data are pro- 
cessed to extract the flux power spectrum. Among these uncertainties the most 
important are: i) mean flux level of the data sets, which could be different from 
data set to data set and is usually kept as a free parameter to marginalize over in 
a reasonable range of observed values (e.g. [127]); ii) the continuum of the distant 
QSOs, which is either subtracted from the data (e.g. [61]) or implicitly removed by 
considering only a smaller portion of the spectra where the continuum is not con- 
tributing to the flux power (e.g. [IOn]); in) the metal absorption systems in the Ly-a 
forest, whose contribution to the flux power could be estimated from Voigt profile 
fitting statistics for high-resolution spectra [108] or modelled using the framework 
proposed by [61] for low-resolution spectra; iv) the contribution to the flux power 
due to the damping wings of strong damped Lj-a systems (DLAs) which again can 
be either removed (VHS) or modelled with the help of hydro simulations (SDSS); v) 
uncertainties in the noise level and resolution of the instrument, which is particularly 
important for the SDSS data set. We stress that among all these uncertainties the 
first two are the most important, while the metal line contribution and the presence 
of DLAs is modelled using parameters (constrained with priors) that are varied in 
the MCMC method. We refer to the extensive analysis made in Refs. [61, 112] for a 
more quantitative discussion on the systematic errors and in the raw data processing 
technique. 

2. Theoretical uncertainties, related to numerical simulations, such as the finite 
box size and finite numerical resolution. These are usually addressed by running a 
number of simulations that explore further the parameter space in resolution and 
box-size and checking the impact on the flux power. Note that these corrections are 
of course code-dependent and could be different in Eulerian and Lagrangian code. 
For a comparison between different codes we refer to [107]. It is also important to 
compare the full hydrodynamical simulations with other less time consuming hydro 
schemes such as the HPM that have been used by several authors. The agreement 
between the Eulerian code ENZO and the Lagrangian code gadget-ii in terms of 
flux power is found to be below or around the 5% level [107], while the discrepancy 
between gadget-ii and HPM [128] could be as high as 20-30% (in the low redshift 
bins, while at higher redshift is usually better and at the percent level) as found 
in [106]. A quantitative analysis of Ly-a flux and line properties using the code 
ENZO has also been made in [129] with a particular focus on box-size and resolution 
problems. 

3. Astrophysical uncertainties relevant for Ly-a physics, e.g. concerning the ther- 
mal history and ionization history of the IGM. One of the key parameters describing 
the influence of the IGM on the flux power spectrum is 7, which appears in the 
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empirical relation T = To(l + S)'^~^ between density fluctuations and the IGM tem- 
perature. Recent results of [■>')] suggested an inverted temperature-density relation 
(i.e. the possibility of having 7 < 1). In order to derive conservative bounds, it was 
suggested in [55] to marginalize 7 over a wide range, including values below 7 = 1. 
In this case, the effect of 7 on the flux power spectrum is not necessarily captured by 
a linear approximation, and one should use at least a second order Taylor expansion 
of Pf(/c, z) as a function of 7 [55, 68]. The results of this paper assume a wide range 
and a second-order Taylor expansion for 7, as in these last references. This issue is 
particularly relevant for AWDM and ACWDM analyses, because the the effect of 
changing the density-temperature relation can be somewhat degenerate with that of 
introducing a cut-off in the primordial matter power spectrum. Indeed, we found 
that implementing the old method used in [65, 67] (7 > 1, flrst order Taylor expan- 
sion of Pp in 7) affects our results for the mass bound m^|p by about 20%. Since the 
density-temperature relation issue is not fully settled yet, we consider the uncertainty 
on 7 to be one important source of systematic errors in our analysis. Another phys- 
ical ingredient which is particularly important and so far has been poorly modelled 
in hydrodynamical simulations is the reionization of hydrogen (which is related to 
the flltering scale of the IGM, see for example [19]): estimates give a few percent im- 
pact on the SDSS flux power but radiative transfer effects (that are likely to generate 
temperature and ionization fluctuations at large scales) need to be better understood 
in order to quantify the effect on the flux power at several scales and redshifts at 
a level comparable to the statistical error of the data [96, 130]. In principle all the 
physical processes that could impact on the flux power at large scales (roughly above 
few Mpc) need to be considered: galactic winds, star formation, cosmic rays. Active 
Galactic Nuclei (AGN) feedback, etc. (see e.g. [109, 60, 112, 131]). However, we 
must stress that: i) accurate investigations made with hydrodynamical simulations 
show that the impact of these effects at large scales is either minimal or degenerate 
with other parameters of the thermal history that are marginalized over; ii) the effect 
of a warm dark matter component on the flux power spectrum does not seem to be 
very degenerate with other physical/cosmological parameters (at least in the AWDM 
case), see e.g. [6G]. 

4. Uncertainties related to the parameter extraction method. As discussed al- 
ready in Sections 5 and 6.1.1, the ambiguity of deflning a stopping criteria for the 
MCMCs (corresponding to a given degree of convergence of the chains) can lead 
to significant errors in the bounds. This error can always be reduced by increas- 
ing the length of the chains, but for models like ACWDM for which the execution 
time of CAMB is larger than usual (due to the necessity of including several discrete 
momenta in the evolution equations of the warm component), reaching a very high 
degree of convergence quickly becomes prohibitive in terms of CPU time. The er- 
ror on the bounds can be estimated by CoSMoMC itself, as already mention in 
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the result sections. It should always be included when quoting results from Cos- 
MoMC. For instance, it lead us to write our AWDM mass bound for the SDSS case 
as m^|p=12.1 keV instead of 13.5 keV. 

5. Finally, there are also uncertainties related to the CMB data. For instance, it 
is well-known that a different treatment of systematic errors in WMAP5 (beam shape, 
small-/ likelihood, etc.) has led to a shift of the allowed range of some parameters 
with respect to WMAP3 (e.g. cxs). However, the results of Fig. 11 show that the 
parameters describing the warm sector are not degenerate with other cosmological 
parameters, and should therefore be rather insensitive to changes in the CMB data. 

In summary, our estimates are based on the current state-of-the art concerning 
the modelling of systematic errors involved in this field. We checked many of these 
assumptions with designated test runs and found that the bounds remain intact 
within ^ 30% uncertainty. Of course, the number of assumptions is such that derived 
bounds should be considered with care. As an illustration of this, the mass bound 
obtained by Seljak et al. [66] and Viel et al. [67] for the same model (pure AWDM) 
and data set (SDSS Lj-a plus CMB data) differ by 30% due to variations in the 
treatment of systematics, data interpretation, numerical issues, etc. 

8 Conclusion 

In this work we have performed a thorough analysis of Lj-a constraints on warm 
and cold plus warm DM models using WMAP5 data, combined with VHS or SDSS 
Ly-a data. As expected from the quality of the data, the results based on SDSS 
Ly-a data are much more constraining, and below we only summarize this case. 

Our results apply to any model in which the phase space distribution of the warm 
component can be approximated by a Fermi-Dirac or rescaled Fermi-Dirac distribu- 
tion. The latter case (assuming a temperature T = T^) is a first-order approximation 
for the case of non-resonantly produced sterile neutrinos. In most of the paper, we 
expressed our mass bounds in terms of the mass jtInrp defined in this case. 

Within a conservative frequentist approach, for pure AWDM models, masses mNpp 
below 8 keV are excluded at 99.7% CL (mNRP > 9.5 keV at 95% CL). In the case of 
CWDM models, our analysis shows that for m-f^^p = 5 keV (the smallest WDM mass 
probed in this investigation) as much as 60% of WDM is allowed at 99.7% CL (40% 
at 95% CL). We believe that these results are robust within the current state-of-the- 
art. Of course, as the Ly-a method is still being developed, these results could be 
affected by still unknown systematic uncertainties (see discussion in the Section 7). 

From the same data we have also obtained Bayesian credible intervals on the same 
parameters. The bound on mNRP (at 95% CL), defined as in previous works [65, 67, 
66]), are m^RP > 12.1 keV for pure AWDM; for ACWDM, the joint probability of 
(-PwDM, l/''^NRp) is displayed in figure 12: this shows that for the smallest probed mass 
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Figure 14: Comparison of X-ray bounds with the results of our Ly-a analysis of 
ACWDM, for the case in which NRP sterile neutrinos contribute to 60% (upper 
panel) or 20% (lower panel) of the dark matter density. The wide gray band (labeled 
"NRP production" ) shows the results of the computation of the NRP sterile neutrino 
abundance, taking into account uncertainties from QCD (see [132, 82] for details). 
The X-ray bounds are from [45] (label "MW"), from [133] (label "M31") and from [41] 
(yellow region in the lower right corner). To account for possible uncertainties in the 
determination of the DM amount in various object, the X-ray bounds have been 
rescaled by a factor of 2 (c.f. [133]). The vertical dashed line marks the 5 keV lower 
bound on the mass found in the present analysis, assuming 60% of WDM. For 20% 
of DM, we found that the mass remains unconstrained by Ly-a data. 

5 keV, the data is compatible with -Fwdm < 35%. 

Our results can be easily translated for the case thermal relics, the corresponding 
limits are the following. In the frequentist framework, for pure AWDM, masses 
below 1.5 keV are excluded at 99.7% CL (mTR > 1.7 keV at 95% CL). In the case 
of CWDM models, our analysis shows that for jtt-tr — 

1.1 keV (the smallest WDM 
mass probed in this investigation) as much as 60% of WDM is allowed at 99.7% CL 
(40% at 95% CL). In the Bayesian approach, for pure AWDM, the 95% CL bound 
defined in the same way as before reads > 2.1 keV. For the ACWDM, the joint 
constraints on the mass and on the warm fraction can be easily obtained from the 
left plot in figure 12, translated in terms of using equation (28). 

Let us now discuss the implications of our results for models in which the dark 
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matter (or part of it) consists of non-resonantly produced sterile neutrinos. 

The lower mass bound for the pure AWDM case should be compared with X-ray 
bounds [43, 44, 41, 134, 135, 4(3, 136, 45, 47, 133], that restrict the mass of NRP 
sterile neutrinos from above. Confronting the upper bound mNRP < 4 keV (see e.g. 
Ref. [133]) with the bound m^^p > 8 keV obtained in this work (both at 99.7% CL) 
excludes the NRP scenario. Again, this conclusion should be taken with care in view 
of the above discussion of systematic errors. 

Next, we analyze the case of mixture of NRP sterile neutrino and some other, cold 
DM particle (similar to [94]). To this end, one should rescale the X-ray constraints 
of [43, 44, 41, 134, 135, 46, 136, 45, 47, 133] and the results of the NRP DM 
abundance computation [132, 82] by a factor F^dm < 1- We can compare the upper 
bounds rescaled in this way with the results of our CWDM runs (see Fig. 14). We 
see that for Fwdm > 0.6, SDSS Lj-a bounds are in conflict with X-ray constraints 
at 99.7% CL. For an NRP fraction 

-fwDM — 0.4 — 0.6, an allowed window for the 
mass appears around ?Tij^pp ~ 5 keV (c.f. Fig. 14, upper panel). For smaller Fwdmi 
the upper mass bound from X-rays quickly increases (as the lower panel of Fig. 14 
demonstrates). At the same time, given that 1 and 2cr Ly-tt contours for -Fwdm 
become horizontal approaches 5 keV (c.f. Fig. 12), we expect that masses 

smaller than 5 keV with Fwdm ^ 0.4 should be allowed (although a precise analysis 
should properly take into account thermal velocities, c.f. Sec. 4). Thus, we conclude 
that for FwDM < 0.6 a window of allowed NRP masses opens up, and for F^dm ~ 0.2 
we see that all masses (above the phase-space density bound [25]) are allowed up to 
~ 15 keV. 

Another DM model which has similarities with the ACWDM scenario discussed 
here is that of resonantly produced sterile neutrinos [80, 85, 83]. In this case, the DM 
velocity dispersion has a colder (resonant) and a warmer (NRP produced) component. 
In this specific analysis is needed in order to apply the results of this work. 

We will present these results elsewhere [137]. 

Further improvements in the results presented here would require performing sim- 
ulations for a larger grid of points (mNRP, -Fwdm), with a better control of systematics 
as well as better understanding of the physics of the IGM. Ideally, one would like to 
fit all the IGM statistics at once (flux power spectrum, flux probability distribution 
function, flux bispectrum and all the line-based statistics) in the spirit of [138, 55] 
to find a consistent IGM model from low to high redshifts. 

High redshift QSOs are important for WDM (and any matter power spectrum 
related) investigations, since at higher redshifts the structures are more linear and 
should more closely resemble the underlying matter power spectrum. However, at z > 
5 the IGM becomes more neutral and the mean flux level lowers. Thus, interpreting 
the almost saturated spectra would require hydrodynamical simulations at much 
higher resolution in order to account for very small structures that contribute to 
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absorption. Furthermore, the number density of QSOs drops significantly at high 
redshifts and thereby the progress will be somewhat limited also in a statistical sense. 

In principle, QSO spectra (especially those at high resolution and with high signal- 
to-noise) carry information down to scales comparable to the resolution element (for 
VLT, Very Large Telescope, UVES-spectra ~ 0.02 A which could be roughly trans- 
lated into 200 comoving kpc at 2 = 3), thereby starting to probe the sub-Mpc 
scales. Although these numbers are very promising, it is important to stress that at 
these scales the flux power spectrum becomes very steep, so that a small change in 
the astrophysical and cosmological parameters could have a large impact on it. To 
interpret it with hydrodynamical simulations could become challenging both for res- 
olution and physical issues (all the physical ingredients that could have an influence 
on galactic scales need to be modelled and their effect re-computed). It is probably 
more promising to tackle some physical effects such as the IGM thermal evolution 
with independent probes (like line statistics) and use them as strong priors in the 
interpretation of the flux power at small scales, in order to restrict the parameter 
space. A large number of intermediate resolution QSO spectra, between that of SDSS 
and VLT spectra, will be also very important in order to cross-check the effects of 
systematics on the data and improve the limits obtained in this work. This effort 
could be possible quite soon with the X-shooter spectrograph [139]. 

Acknowledgments 

We would like to thank M. Laine, A. Maccio, M. Shaposhnikov for help at vari- 
ous stages of this project. We are grateful to V. Savchenko for his invaluable help 
with performing CoSMoMC runs and adapting CoSMoMC to use GRID tech- 
nologies. JL acknowledges support from the EU 6th Framework Marie Curie Re- 
search and Training network "UniverseNet" (MRTN-CT-2006-035863). MV thanks 
ASI/AAE theory grant and INFN-PD51 for financial support. OR acknowledges 
support from the Swiss National Science Foundation. Hydrodynamical simulations 
were performed on the COSMOS supercomputer at DAMTP and at High Perfor- 
mance Computer Cluster Darwin (HPCF) in Cambridge (UK). The Darwin Su- 
percomputer of the University of Cambridge High Performance Computing Service 
(http: / /www. hpc.cam.ac.uk) is provided by Dell Inc. using Strategic Research Infras- 
tructure Funding from the Higher Education Funding Council for England. COS- 
MOS is a UK-CCC facility which is supported by HEFCE, PPARC and Silicon 
Graphics/Cray Research. Most MCMC simulations were performed on the computer 
clusters of the Ukrainian Academic GRID, and in particular on the cluster in the Bo- 
golybov Institute for Theoretical Physics (Kiev) and in the Institute for Scintillation 
Materials (ISMA), Kharkov. We are thankful to E. Martynov and S. Svistunov for 
providing for us the possibility to use these facilities. The remaining part of MCMC 



41 



runs was performed on the MUST cluster at LAPP (IN2P3/CNRS and Universite 
de Savoie). 

Appendices 

A Confidence vs. credible interval 

The main goal of this project is to obtain reliable bounds on the parameters of some 
particle physics DM models, based on the cosmological observations of the Ly-a 
forest. This requires a careful examination of the notion of confidence limits in the 
context of particle physics and cosmology. An extended discussion of this subject 
is beyond the scope of this paper, see e.g. [140], or [141]. For a short introduction 
see [142]. 

In particle physics one can usually repeat an experiment to determine the value 
of a parameter p. The measured data is a fluctuating random variable. Then, for 
example, the 95% confidence interval on a parameter p is defined as the answer to 
the following question: "What is the interval within which the parameter p will lie 
in 95% of all measurements?" This approach (often called the "frequentist" one) is 
well suitable for ruling out models. 

In cosmology one needs to address a different problem, because we can only make 
measurements in our observable Universe, i.e. in a single realization of an underly- 
ing statistical model. In the Bayesian approach one defines, first, the likelihood 
C{D,{pi}) of a data set D given a model with parameters {pi}, and then the Prob- 
ability Density Function (PDF) 11 ({pj}) of the parameter set {pi}, which is equal to 
the likelihood C{D°^^, {Pi}) (seen as a function of the parameters, with the data fixed 
to the observed ones) times the prior distribution of each parameter [143, 144, 145]. 
Thus, ultimately, the parameters (not the data!) are treated as a random variables. 
For each parameter pi, one can construct a probability VsiPi) by marginalizing (i.e., 
integrating) n({pj}) over other parameters Pjj^i, and defines the credible interval as 
the interval of Pi containing 95% of the probability VBiPi)- 

However, it is possible to estimate the frequentist confidence limits, starting from 
the same likelihood function. For each parameter p^, one can maximize the likelihood 
over the other parameters, which gives an (unnormalised) probability distribution 
'PpiPi) = iiiaxp^._^. C{D"^^, {pi}). The 95% confidence limits on pi are defined in such 
way that inside the corresponding interval, VpiPi) is greater than some number such 
that 95% of possible data sets lead to a likelihood greater than this number. 

In general, the two intervals (frequentist and Bayesian) do not coincide. However, 
in cosmology, most authors compute Bayesian intervals and call them confidence 
limit. This is the case of all recent bounds on the WDM mass [65, 67, 66, 68]. 



42 



However, in this paper, we want to carry the Bayesian and frequentist analyses in 
parallel. 

The Monte-Carlo Markov Chain (MCMC) method with Metropolis-Hastings al- 
gorithm [146, 104, 143, 145], implemented in the code CoSMoMC [104], became 
a standard method for exploring the Bayesian PDF H({pj}), which is equal to the 
likelihood C{D°^^, {pi}) (modulo a normalization factor) in case of fiat priors. The 
MCMC method is very effective in scanning H (or C) in the region where this func- 
tion is the largest, with the crucial property that once the chains have converged, 
the density of points in parameter space is proportional to H({pj}). This makes it 
an ideal tool for a Bayesian statistical analysis (see e.g. [1 15] and references therein), 
since the credible interval for a given parameter marginalized over other parameters 
is simply obtained by counting the number of chain points in each parameter bin. 

In order to obtain frequentist confidence limits on a parameter pi, we should take 
into account the exact form of the likelihood as a function of the data (in cosmol- 
ogy, the likelihood is usually non-Gaussian, in particular for small-scale CMB data). 
Instead, we choose to use the same approximation as in many frequentist analy- 
ses: we define a function as x^iPi) = niinp^.__,J— 21n£], search for its minimum 
Xmin = Kii%Jx^(Pj)]) finally define the la, 2a, 3a confidence intervals as the re- 
gion in which Ax^iPi) = (Pi) " Xmm is smaller than 1,4 and 9. In order to minimize 
[— 21n£] over Pjj^i, we still use CoSMoMC with pi kept fixed to different values. It 
is well-known that the CoSMOMC algorithm is not optimal for minimization (with 
respect, e.g., to a simulated annealing method), but in this work it turned out to be 
efficient enough. 

References 

[1] S. Calchi Novati, Microlensing in Galactic Halos, Nuovo Cim. 122B (2007) 
557-567 [0711.4474]. 

[2] Y. B. Zel'dovich, Gravitational instability: An approximate theory for large 
density perturbations., A&A 5 (Mar., 1970) 84-89. 

[3] G. S. Bisnovatyi-Kogan, Cosmology with a nonzero neutrino rest mass, AZh 
57 (Oct., 1980) 899-902. 

[4] J. R. Bond, G. Efstathiou and J. Silk, Massive neutrinos and the large-scale 
structure of the universe, Phys. Rev. Lett. 45 (Dec, 1980) 1980-1984. 

[5] A. G. Doroshkevich, M. I. Khlopov, R. A. Sunyaev, A. S. Szalay and I. B. 
Zeldovich, Cosmological impact of the neutrino rest mass. New York Academy 
Sciences Annals 375 (Dec, 1981) 32-42. 



43 



J. R. Bond and A. S. Szalay, The collisionless damping of density fluctuations 
in an expanding universe, ApJ 274 (Nov., 1983) 443-468. 

L. Bergstrom, Non-baryonic dark matter - observational evidence and 
detection methods, Kept. Prog. Phys. 63 (2000) 793 [hep-ph/0002126]. 

G. Bertone, D. Hooper and J. Silk, Particle dark matter: evidence, candidates 
and constraints, Phys. Rep. 405 (Jan., 2005) 279-390 
[arXiv:hep-ph/0404175]. 

J. Carr, G. Lamanna and J. Lavalle, Indirect detection of dark matter. 
Reports of Progress in Physics 69 (Aug., 2006) 2475-2512. 

M. Taoso, G. Bertone and A. Masiero, Dark Matter Candidates: A Ten-Point 
Test, JCAP 0803 (2008) 022 [0711.4996]. 

S. L. Dubovsky, P. G. Tinyakov and I. I. Tkachev, Massive Graviton as a 
Testable Cold- Dark- Matter Candidate, Phys. Rev. Lett. 94 (May, 2005) 
181102-+ [hep-th/0411158]. 

F. Wilczek, Problem of strong P and T invariance in the presence of 
mstantons, Phys. Rev. Lett. 40 (1978) 279-282. 

S. Weinberg, A new light boson?, Phys. Rev. Lett. 40 (1978) 223-226. 

R. Holman, G. Lazarides and Q. Shafi, Axions and the dark matter of the 
Universe, Phys. Rev. D 27 (Feb., 1983) 995-997. 

S. Dodelson and L. M. Widrow, Sterile-neutrinos as dark matter, Phys. Rev. 
Lett. 72 (1994) 17-20 [hep-ph/9303287]. 

H. Pagels and J. R. Primack, Super symmetry, cosmology, and new physics at 
teraelectronvolt energies, Phys. Rev. Lett. 48 (Jan., 1982) 223-226. 

B. W. Lee and S. Weinberg, Cosmological lower bound on heavy-neutrino 
masses, Phys. Rev. Lett. 39 (1977) 165-168. 

H. E. Haber and G. L. Kane, The search for supersymmetry: Probing physics 
beyond the standard model, Phys. Rep. 117 (Jan., 1985) 75-263. 

L. Covi, J. E. Kim and L. Roszkowski, Axinos as Cold Dark Matter, 
Phys. Rev. Lett. 82 (May, 1999) 4180-4183 [arXiv:hep-ph/9905212]. 

L. Covi, H.-B. Kim, J. E. Kim and L. Roszkowski, Axinos as dark matter, 
JHEP 05 (2001) 033 [hep-ph/0101009]. 



44 



A. Kusenko and M. E. Shaposhnikov, Supersymmetric Q-balls as dark matter, 
Phys. Lett. B418 (1998) 46-54 [hep-ph/9709492]. 

V. A. Kuzmin and I. I. Tkachev, Ultra-high energy cosmic rays, superheavy 
long-living particles, and matter creation after inflation, JETP Lett. 68 
(1998) 271-275 [hep-ph/9802304]. 

D. J. H. Chung, E. W. Kolb and A. Riotto, Superheavy dark matter, 
Phys. Rev. D 59 (Jan., 1999) 023501-+ [arXiv:hep-ph/9802238]. 

S. Tremaine and J. E. Gunn, Dynamical role of light neutral leptons in 
cosmology, Phys. Rev. Lett. 42 (1979) 407-410. 

A. Boyarsky, O. Ruchayskiy and D. lakubovskyi, A lower hound on the mass 
of Dark Matter particles, JCAP, to appear (2008) [0808.3902]. 

J. Madsen and R. I. Epstein, Firm hounds on the neutrino mass from the 
distribution of dark matter in galaxies, ApJ 282 (July, 1984) 11-18. 

J. Madsen, Generalized Tremaine-Gunn limits for bosons and fermions, 
Phys. Rev. D 44 (Aug., 1991) 999-1006. 

J. Madsen, Phase-space constraints on bosonic and fermionic dark matter. 
Physical Review Letters 64 (June, 1990) 2744-2746. 

J. J. Dalcanton and C. J. Hogan, Halo cores and phase space densities: 
Observational constraints on dark matter physics and structure formation, 
ApJ 561 (2001) 35-45 [astro-ph/0004381]. 

C. J. Hogan and J. J. Dalcanton, New dark matter physics: Glues from halo 
structure, Phys. Rev. D62 (2000) 063511 [astro-ph/0002330]. 

J. Madsen, Dark matter phase space densities, Phys. Rev. D 64 (July, 2001) 
027301-+ [arXiv : astro-ph/0006074] . 

J. Madsen, Dark matter phase space densities, Phys. Rev. D 64 (July, 2001) 
027301-+ [arXiv : astro-ph/0006074] . 

D. Boyanovsky, H. J. de Vega and N. G. Sanchez, Gonstraints on dark matter 
particles from theory, galaxy observations, and N-body simulations, 

Phys. Rev. D 77 (Feb., 2008) 043518-+ [arXiv: 0710 . 5180]. 

D. Gorbunov, A. Khmelnitsky and V. Rubakov, Gonstraining sterile neutrino 
dark matter by phase-space density observations, JGAP 0810 (2008) 041 
[0808.3910]. 



45 



E. Kolb and M. Turner, The Early Universe. Addison- Wesley, Reading, MA, 
USA, 1990. Prepared with KTeX. 

J. L. Feng, A. Rajaraman and F. Takayama, Superweakly-interacting massive 
particles, Phys. Rev. Lett. 91 (2003) 011302 [hep-ph/0302215]. 

F. Takayama and M. Yamaguchi, Gravitino dark matter without R-parity, 
Phys. Lett. B485 (2000) 388-392 [hep-ph/0005214]. 

W. BuchmuUer, L. Covi, K. Hamaguchi, A. Ibarra and T. Yanagida, 
Gravitino dark matter in R-parity breaking vacua, JHEP 03 (2007) 037 
[hep-ph/0702184]. 

J. P. Conlon and F. Quevedo, Astrophysical and Gosmological Implications of 
Large Volume String Gompactifications, JGAP 0708 (2007) 019 [0705.3460]. 

M. Lattanzi and J. W. F. Valle, Decaying warm dark matter and neutrino 
masses, Phys. Rev. Lett. 99 (2007) 121301 [0705.2406]. 

A. Boyarsky, A. Neronov, O. Ruchayskiy, M. Shaposhnikov and I. Tkachev, 
How to find a dark matter sterile neutrino?, Phys. Rev. Lett. 97 (dec, 2006) 
261302 [astro-ph/0603660]. 

G. Bertone, W. Buchmuller, L. Covi and A. Ibarra, Gamma-rays from 
decaying dark matter, JGAP 0711 (2007) 003 [0709.2299]. 

A. Boyarsky, A. Neronov, O. Ruchayskiy and M. Shaposhnikov, Gonstraints 
on sterile neutrino as a dark matter candidate from the diffuse X-ray 
background, MNRAS 370 (July, 2006) 213-218 [astro-ph/0512509]. 

A. Boyarsky, A. Neronov, O. Ruchayskiy and M. Shaposhnikov, Restrictions 
on parameters of sterile neutrino dark matter from observations of galaxy 
clusters, Phys. Rev. D 74 (nov, 2006) 103506 [astro-ph/0603368]. 

A. Boyarsky, J. Nevalainen and O. Ruchayskiy, Gonstraints on the parameters 
of radiatively decaying dark matter from the dark matter halo of the milky 
way and ursa minor, A&A 471 (Aug., 2007) 51-57 [astro-ph/0610961]. 

A. Boyarsky, O. Ruchayskiy and M. Markevitch, Gonstraints on parameters 
of radiatively decaying dark matter from the galaxy cluster le0657-56, ApJ 
673 (2008) 752 [astro-ph/0611168]. 

A. Boyarsky, J. W. den Herder, A. Neronov and O. Ruchayskiy, Search for 
the light dark matter with an x-ray spectrometer, Astropart. Phys. 28 (2007) 
303-311 [astro-ph/0612219]. 



46 



[48] P. Bode, J. P. Ostriker and N. Turok, Halo formation in warm dark matter 
models, ApJ 556 (2001) 93-107 [astro-ph/0010389]. 

[49] L. Hui, N. Y. Gnedin and Y. Zhang, The Statistics of Density Peaks and the 
Column Density Distribution of the Ly alpha Forest, ApJ 486 (Sept., 1997) 
599-+ [astro-ph/9608157]. 

[50] N. Y. Gnedin and A. J. S. Hamilton, Matter power spectrum from the 

lyman- alpha forest: Myth or reality?, Mon.Not.Roy.Astron.Soc. 334 (2002) 
107-116 [astro-ph/0111194]. 

[51] D. H. Weinberg, R. Dave, N. Katz and J. A. Kollmeier, The Lyman-a Forest 
as a Cosmological Tool, in AIP Conf. Proc. 666: The Emergence of Cosmic 
Structure (S. H. Holt and C. S. Reynolds, eds.), pp. 157-169, May, 2003. 

[52] A. A. Meiksin, The Physics of the Intergalactic Medium, ArXiv e-prints 
(Nov., 2007) [0711.3358]. 

[53] M. Viel, M. G. Haehnelt, R. F. Carswell and T. S. Kim, The effect of (strong) 
discrete absorption systems on the lyman alpha forest flux power spectrum, 
Mon. Not. Roy. Astron. Soc. 349 (2004) L33 [astro-ph/0308078]. 

[54] T. S. Kim, J. S. Bolton, M. Viel, M. G. Haehnelt and R. F. Carswell, An 
improved measurement of the flux distribution of the Ly-alpha forest in QSO 
absorption spectra: the effect of continuum fitting, metal contamination and 
noise properties, 0711.1862. 

[55] J. S. Bohon, M. Viel, T. S. Kim, M. G. Haehnelt and R. F. Carswell, Possible 
evidence for an inverted temperature- density relation in the intergalactic 
medium from the flux distribution of the Lyman-alpha forest, 0711 . 2064. 

[56] L. Gao and T. Theuns, Lighting the Universe with Filaments, Science 317 
(Sept., 2007) 1527- [arXiv: 0709. 2165]. 

[57] P. L. Biermann and A. Kusenko, Relic keV sterile neutrinos and reionization, 
Phys. Rev. Lett. 96 (2006) 091301 [astro-ph/0601004]. 

[58] J. Stasielak, P. L. Biermann and A. Kusenko, Thermal evolution of the 

primordial clouds in warm dark matter models with kev sterile neutrinos, ApJ 
654 (Jan., 2007) 290-303 [arXiv:astro-ph/0606435]. 

[59] C. . Faucher-Giguere, J. X. Prochaska, A. Lidz, L. Hernquist and 
M. Zaldarriaga, A Direct Precision Measurement of the Intergalactic 
Lyman-alpha Opacity at 2<2<4.2, 0709.2382. 



47 



[60] M. Viel, M. G. Haehnelt and V. Springel, Inferring the dark matter power 
spectrum from the Lyman- alpha forest in high-resolution QSO absorption 
spectra, Mon. Not. Roy. Astron. Soc. 354 (2004) 684 [astro-ph/0404600]. 

[61] P. McDonald, U. Seljak, S. Buries, D. J. Schlegel, D. H. Weinberg, R. Gen, 
D. Sliih, J. Schaye, D. P. Schneider, N. A. Bahcall, J. W. Briggs, 
J. Brinkmann, R. J. Brunner, M. Fukugita, J. E. Gunn, Z. Ivezic, S. Kent, 
R. H. Lupton and D. E. Vanden Berk, The Lya Forest Power Spectrum from 
the Sloan Digital Sky Survey, ApJS 163 (Mar., 2006) 80-109 
[arXiv:astro-ph/0405013]. 

[62] M. Viel and M. G. Haehnelt, Cosmological and astrophysical parameters from 
the SDSS flux power spectrum and hydrodynamical simulations of the 
Lyman- alpha forest, Mon. Not. Roy. Astron. Soc. 365 (2006) 231-244 
[astro-ph/0508177]. 

[63] S. H. Hansen, J. Lesgourgues, S. Pastor and J. Silk, Closing the window on 
warm dark matter, Mon. Not. Roy. Astron. Soc. 333, 544 (2002) 
[astro-ph/0106108]. 

[64] K. Abazajian, Linear cosmological structure limits on warm dark matter, 
Phys. Rev. D71 73 (2006) 063513 [astro-ph/0512631] 

[65] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese and A. Riotto, 

Constraining warm dark matter candidates including sterile neutrinos and 
light gravitinos with wmap and the Lyman- alpha forest, Phys. Rev. D71 
(2005) 063534 [astro-ph/0501562]. 

[66] U. Seljak, A. Makarov, P. McDonald and H. Trac, Can sterile neutrinos he 
the dark matter?, Phys. Rev. Lett. 97 (2006) 191303 [astro-ph/0602430]. 

[67] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese and A. Riotto, Can 
sterile neutrinos he ruled out as warm dark matter candidates?, Phys. Rev. 
Lett. 97 (2006) 071301 [astro-ph/0605706]. 

[68] M. Viel, G. D. Becker, J. S. Bolton, M. G. Haehnelt, M. Ranch and W. L. W. 
Sargent, How cold is cold dark matter? Small scales constraints from the flux 
power spectrum of the high-redshift Lyman-alpha forest, Phys. Rev. Lett. 100 
(2008) 041304 [0709.0131]. 

[69] J. L. Feng, S.-f. Su and F. Takayama, SuperWIMP gravitino dark matter from 
slepton and sneutrino decays, Phys. Rev. D70 (2004) 063514 
[hep-ph/0404198]. 



48 



[70] L. Roszkowski, R. Ruiz de Austri and K.-Y. Choi, Gravitino dark matter in 
the CMSSM and implications for leptogenesis and the LHC, JHEP 08 (2005) 
080 [hep-ph/0408227]. 

[71] D. G. Cerdeno, K.-Y. Choi, K. Jedamzik, L. Roszkowski and R. Ruiz de 

Austri, Gravitino dark matter in the GMSSM with improved constraints from 
BBN, JGAP 0606 (2006) 005 [hep-ph/0509275]. 

[72] K.-Y. Choi, J. E. Kim, H. M. Lee and O. Seto, Neutralino dark matter from 
heavy axmo decay, Phys. Rev. D77 (2008) 123501 [0801.0491]. 

[73] M. Bolz, A. Brandenburg and W. BuchmuUer, Thermal production of 
gravitinos, Nucl.Phys. B 606 (2001) 518-544 [hep-ph/0012052]. 

[74] J. Pradler and F. D. Steffen, Thermal Gravitino Production and Gollider 
Tests of Leptogenesis, Phys. Rev. D75 (2007) 023509 [hep-ph/0608344]. 

[75] V. S. Rychkov and A. Strumia, Thermal production of gravitinos, Phys. Rev. 
D75 (2007) 075011 [hep-ph/0701104]. 

[76] S. Borgani, A. Masiero and M. Yamaguchi, Light gravitinos as mixed dark 
matter, Phys. Lett. B386 (1996) 189-197 [hep-ph/9605222]. 

[77] T. Asaka, S. Blanchet and M. Shaposhnikov, The numsm, dark matter and 
neutrino masses, Phys. Lett. B631 (2005) 151-156 [hep-ph/0503065]. 

[78] T. Asaka and M. Shaposhnikov, The nuMSM, dark matter and baryon 

asymmetry of the universe [rapid communication], Phys. Lett. B 620 (July, 
2005) 17-26 [arXiv:hep-ph/0505013]. 

[79] M. Shaposhnikov, Is there a new physics between electroweak and Planck 
scales?, 0708.3550. 

[80] X.-d. Shi and G. M. Fuller, A new dark matter candidate: Non-thermal sterile 
neutrinos, Phys. Rev. Lett. 82 (1999) 2832-2835 [astro-ph/9810076]. 

[81] K. Abazajian, G. M. Fuller and M. Patel, Sterile neutrino hot, warm, and 
cold dark matter, Phys. Rev. D 64 (2001) 023501 [astro-ph/0101524]. 

[82] T. Asaka, M. Laine and M. Shaposhnikov, Lightest sterile neutrino 
abundance within the numsm, JHEP 01 (2007) 091 [hep-ph/0612182]. 

[83] M. Laine and M. Shaposhnikov, Sterile neutrino dark matter as a 

consequence of vMSM-induced lepton asymmetry, JGAP Q (June, 2008) 31 — h 
[arXiv: 0804. 4543]. 



49 



[84] A. D. Dolgov and S. H. Hansen, Massive sterile neutrinos as warm dark 
matter, Astropart. Phys. 16 (2002) 339-344 [hep-ph/0009083]. 

[85] M. Shaposhnikov, The nuMSM, leptonic asymmetries, and properties of 
singlet fermions, JHEP 08 (2008) 008 [0804.4542]. 

[86] A. Boyarsky, O. Ruchayskiy and M. Shaposhnikov, The role of sterile 
neutrinos in cosmology and astrophysics, 0901.0011. 

[87] M. Shaposhnikov and I. Tkachev, The numsm, inflation, and dark matter, 
Phys. Lett. B639 (2006) 414-417 [hep-ph/0604236]. 

[88] A. Kusenko, Sterile neutrinos, dark matter, and the pulsar velocities in models 
with a higgs singlet, Phys. Rev. Lett. 97 (2006) 241301 [hep-ph/0609081]. 

[89] K. Petraki and A. Kusenko, Dark-matter sterile neutrinos in models with a 
gauge singlet in the Higgs sector, Phys. Rev. D77 (2008) 065014 [0711.4646]. 

[90] K. Petraki, Small-scale structure formation properties of chilled sterile 
neutrinos as dark matter, Phys. Rev. D 77 (May, 2008) 105004 — h 
[arXiv: 0801. 3470]. 

[91] D. Boyanovsky, Clustering properties of a sterile neutrino dark matter 
candidate, Phys. Rev. D78 (2008) 103505 [0807.0646]. 

[92] S. Khahl and O. Seto, Sterile neutrino dark matter in B — L extension of the 
standard model and galactic 511 keV line, JCAP 0810 (2008) 024 
[0804.0336]. 

[93] J. Wu, C. M. Ho and D. Boyanovsky, Sterile neutrinos produced near the EW 
scale /.■ mixing angles, MSW resonances and production rates, [0902.4278]. 

[94] A. Palazzo, D. Cumberbatch, A. Slosar and J. Silk, Sterile neutrinos as 
subdominant warm dark matter, arXiv : 0707 . 1495 [astro-ph] . 

[95] M. Viel, M. G. Haehnelt and V. Springel, Inferring the dark matter power 
spectrum from the Lyman- alpha forest in high-resolution QSO absorption 
spectra, Mon. Not. Roy. Astron. Soc. 354 (2004) 684 [astro-ph/0404600]. 

[96] P. McDonald et. ai. The linear theory power spectrum from the lyman- alpha 
forest in the sloan digital sky survey, Astrophys. J. 635 (2005) 761-783 
[astro-ph/0407377]. 

[97] S. Colombi, S. Dodelson and L. M. Widrow, Large scale structure tests of 
warm dark matter, Astrophys. J. 458 (1996) 1 [astro-ph/9505029]. 



50 



[98] J. Lesgourgues and S. Pastor, Massive neutrinos and cosmology, Phys. Kept. 
429 (2006) 307-379 [astro-ph/0603494]. 

[99] A. Lewis, A. Challinor and A. Lasenby, Efficient computation of CMB 
anisotropics in closed FRW models, Astrophys. J. 538 (2000) 473-476 
[astro-ph/9911177]. 

[100] W. Hu, D. J. Eisenstein and M. Tegmark, Weighing neutrinos with galaxy 
surveys, Phys. Rev. Lett. 80 (1998) 5255-5258 [astro-ph/9712057]. 

[101] H. Bi Astrophys. J. 405 (1993) 479. 

[102] M. Viel, S. Matarrese, H. J. Mo, M. G. Haehnelt and T. Theuns, Probing the 
Intergalactic Medium with the Lyman alpha forest along multiple lines of sight 
to distant QSOs, Mon. Not. Roy. Astron. Soc. 329 (2002) 848 
[astro-ph/0105233]. 

[103] M. Zaldarriaga, R. Scoccimarro and L. Hui, Inferring the Linear Power 
Spectrum from the Lyman-alpha Forest, Astrophys. J. 590 (2003) 1-7 
[astro-ph/0111230]. 

[104] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: a 
Monte- Carlo approach, Phys. Rev. D66 (2002) 103511 [astro-ph/0205436]. 

[105] T. Theuns, A. Leonard, G. Efstathiou, F. R. Pearce and P. A. Thomas, 
P3M-SPH simulations of the Lyman-alpha Forest, Mon. Not. Roy. Astron. 
Soc. 301 (1998) 478-502 [astro-ph/9805119]. 

[106] M. Viel, M. G. Haehnelt and V. Springel, Testing the accuracy of the 

hydro-pm approximation in numerical simulations of the lyman-alpha forest, 
Mon. Not. Roy. Astron. Soc. 367 (2006) 1655-1665 [astro-ph/0504641]. 

[107] J. A. Regan, M. G. Haehnelt and M. Viel, Numerical Simulations of the 
Lyman-alpha forest - A comparison of Cadget-2 and Enzo, Mon. Not. Roy. 
Astron. Soc. 374 (2007) 196-205 [astro-ph/0606638]. 

[108] T. S. Kim, M. Viel, M. G. Haehnelt, R. F. Carswell and S. Cristiani, The 

power spectrum of the flux distribution in the Lyman- alpha forest of a Large 
sample of UVES QSO Absorption Spectra (LUQAS), Mon. Not. Roy. Astron. 
Soc. 347 (2004) 355 [astro-ph/0308103]. 

[109] R. A. C. Croft et. ai. Towards a precise measurement of matter clustering: 
Lyman- alpha forest data at redshifts 2-4, Astrophys. J. 581 (2002) 20-52 
[astro-ph/0012324]. 



51 



[110] N. Y. Gnedin and L. Hui, Probing the Universe with the Lyman- alpha Forest: 
I. Hydrodynamics of the Low Density IGM, Mon. Not. Roy. Astron. Soc. 296 
(1998) 44-55 [astro-ph/9706219]. 

[Ill] M. Bernardi et al. [SDSS Collaboration], A feature at z 3.2 in the evolution 
of the Ly-alpha forest optical depth, Astron. J. 125 (2003) 32 
[arXiv : astro-ph/0206293] . 

[112] P. McDonald, U. Seljak, R. Cen, P. Bode and J. P. Ostriker, Physical effects 
on the Lya forest flux power spectrum: damping wings, ionizing radiation 
fluctuations and galactic winds, MNRAS 360 (July, 2005) 1471-1482 
[arXiv : astro-ph/0407378] . 

[113] E. Bertschinger, COSMICS: Cosmological Initial Conditions and Microwave 
Anisotropy Codes, astro-ph/9506070. 

[114] A. Klypin and J. Holtzman, Particle- Mesh code for cosmological simulations, 
astro-ph/9712217. 

[115] A. Klypin, Numerical simulations in cosmology i: Methods, 
astro-ph/0005502. 

[116] V. Springel, The cosmological simulation code CADGET-2, Mon. Not. Roy. 
Astron. Soc. 364 (2005) 1105-1134 [astro-ph/0505010]. 

[117] P. Colin, O. Valenzuela and V. Avila- Reese, On the Structure of Dark Matter 
Halos at the Damping Scale of the Power Spectrum with and without Relict 
Velocities, Astrophys. J. 673 (2008) 203-214 [0709.4027]. 

[118] WMAP Collaboration, J. Dunkley et. al. Five- Year Wilkinson Microwave 
Anisotropy Probe (WMAP) Observations: Likelihoods and Parameters from 
the WMAP data, 0803.0586. 

[119] ACBAR Collaboration, C.-l. Kuo et. al.. High resolution observations of the 
cmb power spectrum with acbar, Astrophys. J. 600 (2004) 32-51 
[astro-ph/0212289]. 

[120] A. C. S. Readhead et. al. Extended mosaic observations with the cosmic 

background imager, Astrophys. J. 609 (2004) 498-512 [astro-ph/0402359]. 

[121] F. Piacentini et. al., A measurement of the polarization-temperature angular 
cross power spectrum of the cosmic microwave background from the 2003 
flight of boomerang, ApJ 647 (Aug., 2006) 833-839 [astro-ph/0507507]. 



52 



W. C. Jones et. al, A measurement of the angular power spectrum of the cmb 
temperature anisotropy from the 2003 flight of boomerang, ApJ 647 (Aug., 
2006) 823-832 [astro-ph/0507494]. 

T. E. Montroy et. ai, A measurement of the cmb spectrum from the 2003 
flight of boomerang, ApJ 647 (Aug., 2006) 813-822 [astro-ph/0507514]. 

C. J. MacTavish et. ah, Cosmological parameters from the 2003 flight of 
boomerang, Astrophys. J. 647 (2006) 799 [astro-ph/0507503]. 

M. Tegmark et. ai, Cosmological constraints from the sdss luminous red 
galaxies, astro-ph/0608632. 

P. Astier et. ai. The supernova legacy survey: Measurement of ujm, ^^^a O'nd w 
from the first year data set, Astron. Astrophys. 447 (2006) 31-48 
[astro-ph/0510447]. 

U. Seljak, P. McDonald and A. Makarov, Cosmological constraints from the 
cosmic microwave background and Lyman a forest revisited, MNRAS 342 
(July, 2003) L79-L84 [arXiv:astro-ph/030257l]. 

N. Y. Gnedin and L. Hui, Probing the Universe with the Lyalpha forest - I. 
Hydrodynamics of the low-density intergalactic medium, MNRAS 296 (May, 
1998) 44-55 [arXiv:astro-ph/9706219]. 

D. Tytler, P. Paschos, D. Kirkman, M. L. Norman and T. Jena, How 
simulated spectra of the lya forest change with the size of the box of numerical 
simulations, 0711 . 2529. 

R. A. C. Croft, Ionizing Radiation Fluctuations and Large-Scale Structure in 
the Lya Forest, ApJ 610 (Aug., 2004) 642-662 [arXiv:astro-ph/0310890]. 

M. Jubelgas, V. Springel, T. EnBlin and C. Pfrommer, Cosmic ray feedback in 
hydrodynamical simulations of galaxy formation, A&A 481 (Apr., 2008) 
33-63. 

T. Asaka, M. Laine and M. Shaposhnikov, On the hadronic contribution to 
sterile neutrino production, JHEP 06 (2006) 053 [hep-ph/0605209]. 

A. Boyarsky, D. lakubovskyi, O. Ruchayskiy and V. Savchenko, Constraints 
on decaying dark matter from XMM-Newton observations of M31., MNRAS 
387 (2008) 1361-1373 [arXiv : 0709 . 2301]. 



53 



[134] S. Riemer-S0rensen, S. H. Hansen and K. Pedersen, Sterile Neutrinos in the 
Milky Way: Observational Constraints, ApJ 644 (June, 2006) L33-L36 
[astro-ph/0603661]. 

[135] C. R. Watson, J. F. Beacom, H. Yuksel and T. P. Walker, Direct x-ray 

constraints on sterile neutrino warm dark matter, Phys. Rev. D74 (2006) 
033009 [astro-ph/0605424]. 

[136] K. N. Abazajian, M. Markevitch, S. M. Koushiappas and R. C. Hickox, 
Limits on the Radiative Decay of Sterile Neutrino Dark Matter from the 
Unresolved Cosmic and Soft X-ray Backgrounds, Phys. Rev. D 75 (Mar., 
2007) 063511-+ [arXiv:astro-ph/0611144]. 

[137] A. Boyarsky, J. Lesgourgues, O. Ruchayskiy and M. Viel, Realistic sterile 

neutrino dark matter with keV mass does not contradict cosmological hounds, 
0812.3256. 

[138] D. Tytler, D. Kirkman, J. M. O'Meara, N. Suzuki, A. Orin, D. Lubin, 

P. Paschos, T. Jena, W.-C. Lin, M. L. Norman and A. Meiksin, Cosmological 
parameters a^, the baryon density Ub, the Vacuum Energy Density uoiy, the 
Hubble Constant and the UV Background Intensity from a Calibrated 
Measurement of H I Lya absorption at z=1.9, ApJ 617 (Dec, 2004) 1-28 
[arXiv : astro-ph/0403688] . 

[139] L. Kaper, S. D'Odorico, F. Hammer, R. Pallavicini and P. Kjaergaard 

Rasmussen, X-shooter: a medium-resolution, wide-band spectrograph for the 
VLT, ArXiv e-prints (Mar., 2008) [0803.0609]. 

[140] D. J. C. MacKay, Information Theory, Inference and Learning Algorithms. 
Cambrdige University Press, 2003. 

http : //www. inference .phy . cam. ac .uk/mackay/ itprnn/book.html. 

[141] G. J. Feldman and R. D. Cousins, A Unified approach to the classical 
statistical analysis of small signals, Phys. Rev. D57 (1998) 3873-3889 
[physics/9711021]. 

[142] D. Karlen, Credibility of confidence intervals, . Prepared for Conference on 
Advanced Statistical Techniques in Particle Physics, Durham, England, 18-22 
Mar 2002. 

[143] A. Kosowsky, M. Milosavljevic and R. Jimenez, Efficient cosmological 

parameter estimation from microwave background anisotropics, Phys. Rev. 
D66 (2002) 063007 [astro-ph/0206014]. 



54 



[144] J. Dunkley, M. Bucher, P. G. Ferreira, K. Moodley and C. Skordis, Fast and 
reliable mcmc for cosmological parameter estimation, Mon. Not. Roy. Astron. 
Soc. 356 (2005) 925-936 [astro-ph/0405462]. 

[145] R. Trotta, Bayes in the sky: Bayesian inference and model selection in 
cosmology, Contemporary Phys. 49 (2008) 71-104 [0803.4089]. 

[146] R. M. Neal, "Probabilistic inference using Markov Chain Monte Carlo 
methods." http://cosmologist.info/Neal93, 1993. 



55 



